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1. Introduction. The class of vector differential equations considered in this paper 
is represented by 
d'r/dt’ — f(r) = b(0), (1) 


in which r is a radius vector in a Cartesian frame, f(r) is a point function derivable as the 
gradient of a scalar, and b(¢) is an arbitrary function of the independent variable ¢. In 
its physical application, Eq. (1) determines the position of a particle subject to the 
acceleration b + f at time ¢. To allow the wider interpretation of Eq. (1) as a curve in 
three-space (for fixed initial conditions), ¢ will not be identified with time (except in the 
applications of Sec. 3). The basic assumption will be made that f(r) (which, in general, 
is a non-linear function of components of r) is a slowly varying function as compared 
to are length on any integral curve of Eq. (1). The scalar corresponding to f(r) will be 
assumed differentiable to fourth order at least, and the function b(¢) will be assumed 
differentiable at the initial point of any integral curve of Eq. (1) to third order at least 
(but otherwise merely Riemann-integrable). 

In a previous paper (referred to hereafter as J) by the author and others [1], a special 
case of Eq. (1) was treated, in which f(r) and b(t) were the gravitational and non-gravi- 
tational accelerations respectively of a particle. A quasi-linear approximation for two- 
dimentional motion was described, which takes into account terms in f of higher order 
than does a linearized approximation, and which depends on the assumption that the 
radius of curvature of the corresponding trajectory be slowly varying. In this paper, 
the quasi-linear approximation in question is generalized in two directions: (a) to three 
dimensions (by assuming that the radius of torsion of an integral curve is likewise slowly 
varying); (b) to include an arbitrary f(r) derivable as the gradient of a non-harmonic 
scalar. The approximation exploits the fact that coordinate distances (canonical co- 
ordinates) along the tangent, principal normal, and binormal of an analytic integral 
curve are proportional (within dominant terms) to the first, second, and third powers 
of the are length, respectively, by virtue of the Frenet-Serret formulas. Hence, non- 
linear terms in f through third order in the are length can be expressed as linear functions 
of the canonical coordinates and thus of any Cartesian coordinates. In a sense, the 
approximation considered is the analogue in the non-linear case of Liouville’s approxi- 
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Wentzel-Kramers-Brillouin-Jeffries approximation 


mation [2, 3] (referred to as the 
approximation 


[4, 5, 6] in quantum mechanics). A prefatory discussion of the linearized 
to Eq. (1) is given first. 

Some remarks on the practical use of approximations of the type considered are 
in order. The general function of such approximations is to determine corrections to a 
zero-or first-order solution for a limited range of the independent variable. Approxima- 
tions of this type cannot compete, of course, with modern numerical methods of solution 
for an extended range of the independent variable. For moderate ranges of the inde- 
pendent variable, however, linear approximations can be applied piece-wise by identi- 
fying initial conditions on one interval with terminal conditions on a preceding interval. 
Such a continuation solution has been used in J for the computation of an actual tra- 
jectory. 

2. Linearized Approximation. The scalar potential g, in terms of which f = Veg, 
will be expressed as a function of three coordinates x; (¢ = 1,2,3) forming a Cartesian 
frame (right-handed in the order given) with origin at the initial point O of an integral 
curve in question. The Taylor expansion of ¢ in these coordinates at O will have the 
form 

g=o ° +e7° +e7° +e? 4+°::, (2) 
0 at O and each term is a rational integral homogeneous function of the 2, 


where ¢ = 
yg will be 


(the superscript notation indicates the degree). Through quadratic terms, 
written 

g= DL Await DB 22;, (i,j = 1, 2,3), (3) 
where B;; = B;, . 

' { 


A matrix notation {a} will be used interchangeably with a for a vector, where {a} 


stands for the column vector {a,, a , a3} of components. The linearized approximation 


f‘” to f consists in retaining only terms in ¢ through second order. From Eq. (3), it is 
{ff } =A+ Pir}, (4) 
where {r} = {z,}, A = {A;}, and the coefficient matrix P is 
P = 2B (5) 


in terms of the matrix B 


= [B,,]. With the notation r° = dr/dt, the corresponding linear- 
ized approximation to Eq. (1) is 
ir} — Pir} = {bd} +A, (6) 

which assumes that f is a slowly-varying function of position with respect to are length 
on an integral curve. 

Since P is a symmetric matrix, its eigenvalues are real, and Eq. (6) can be solved 
directly by reduction to normal coordinates. The solution for the normal coordinate 
x’ corresponding to the negative eigenvalue — w; of P is 


vt 


zi = vw, sinw,;t + w; | [b'(7) + Af] sin w,(t — 7) dr, 


where v! is the initial value of dx’/dt and b’ + A’ is the transform in the normal coordi- 


nates of b + A. For a positive eigenvalue of P, the circular functions in Eq. (7) must 
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be replaced by the corresponding hyperbolic functions. In the classical theory of small 
oscillations of a particle, eigenvalues leading to such hyperbolic solutions are usually 
excluded. Such exclusion is not necessary (at least on purely mathematical grounds) 
in the application to particle motion of the linear approximations discussed here, since 
the application is to the initial phase of a forced motion. For this reason, the approxima- 
tions do not necessarily carry implications on the boundedness in the large of the corre- 
sponding motion. Note that, from Earnshaw’s theorem [7], all the eigenvalues of P can 
be of one sign only if g is non-harmonic. 

Ordinarily, the coefficients B;; must be determined by direct expansion of yg. The 
Appendix describes a method of determining these coefficients which relates them to 
geometric parameters of the equipotential passing through the origin O. For the purpose 
of formulating this result, take the x,-axis in the direction of — Vg, and take the 2,—, 
%,— axes in the directions of the lines of curvature of the equipotential surface ¢ = 0 
in the tangent plane of this equipotential at O. Let r, and r, be the principal radii of 
normal curvature of this surface at O in the directions of x, and zx, respectively’. Let 
r,.. be the derivative evaluated at O of the principal radius of normal curvature in the 
x-direction with respect to are length on the line of curvature in the z,-direction; and 


let 72,, be defined correspondingly. The matrix B is then 


a 
r; 6 T2172 

J —§ 1 *] 
B=‘ 0 To 1; 07 ; (8) 

-1 l 1 1 ° 

=f, 7, —T 1.09} —(r7, +r. t+ Go/fo) 

where oy is a mass density defined by Poisson’s equation V*¢ = —o at O, and fy = | Ve | 
at O. The vector A is {0, 0, — fo} in the particular frame introduced here. Equation 


(8) is convenient when the form of the equipotential is a datum given directly (as it is 
for f a gravitational acceleration). Use will be made of it in an application in Sec. 3. 

3. Quasi-Linear Approximation. [or a non-vanishing initial dr/dt, the corresponding 
integral curve is analytic at the initial point O, since V¢ and b are assumed differentiable 
at O (to third order, at least). Hence, the radius vector r to a point on the curve can be 


expanded in the Taylor series [8] 


8 A ( 'B 3° 
r=as+3 o- 148g tie. .-., (9) 
<p p p pt) 6 


where s is the are length on the curve; e, 8, y are the unit tangent, unit principal normal, 
and unit binormal at O respectively; and p, 7, p’ are the (non-vanishing) radius of curva- 
ture, (non-vanishing) radius of torsion, and derivative of p with respect to s respectively 
on the integral curve at O. Equation (9) is equivalent to the standard Frenet-Serret 
formulas. It is understood that a, 8, y, p, 7, p’ are constants evaluated at O, but speci- 
fying subscripts will be omitted for convenience. 

The coefficients of the terms of order s* and higher in Eq. (9) contain products of 
powers and derivatives with respect to s of 1/p and 1/r. If all derivatives of p and + 
with respect to s are sufficiently small, one can take s sufficiently small so that terms 


'The principal radii r; and rz are taken positive or negative according as the corresponding centers 


of curvature lie on the positive or negative half of the 23-axis. 
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of the form s*/p*, s*/p’r, --- (which do not involve derivatives) are negligible in Eq. 
(9). In this case, r is given by the terms of the series (9) which are indicated explicitly. 
The permissible range of s for this approximation is given by s < p, 7 and is thus smaller 
as p and 7 are smaller. Correspondingly, the terms of order s° and s* in Eq. (9) are then 
smaller for p and 7 smaller. Since these terms form the basis of the quasi-linear approxi- 
mation to be described, the range of s over which this approximation yields a significant 
correction to the linearized approximation is greater for p and r larger. 
On these assumptions, Eq. (9) determines any 2; as a linear function of s, s° and s°, 
which is given by 
ir} = F’S{s'}, (10) 


where {s*} designates the column vector, {s, s°, s°}, I’ is the transpose of the matrix I 
of direction cosines of a, 8, y, 


a Qe Qs 


and the matrix S is 


] 0 — (6p) 
B= 10 (2p)~* —p'(6p)' |}. (12) 
0 0 —(6pr)' 


Note that the components of the column vector S {s'} in Eq. (10) are simply canonical 


coordinates of the curve. Equation (10) can be inverted and written 


{s'} = ST {r}, (13) 
where 
l 0 —tp 
s'=1,0 2p —2p’r |, (14) 
0 0 —6pr 


provided S is non-singular, which is true if the determinant 

| S| = —(12p*7)~’ (15) 
be non-vanishing. Thus the following results apply to an integral curve of Eq. (1) which 
is curved (1/p # 0) and twisted (1/7 # 0) at the origin; the special case of an integral 
curve which is plane at the origin will be treated later. Under these restrictions, Eq. 
(13) expresses s, s° and s° linearly in terms of coordinates zx, of the curve. 


If quintic and higher terms in ¢ are neglected, f can be written 


{ff} =A+ Pir} + {Vo} + [Ve }. (16) 
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In general, the cubic and quartic terms of g, yg’ and y™’, are polynomials having re- 
spectively 10 and 15 coefficients (these two sets of coefficients can depend on only 7 
and 9 independent constants respectively if g is harmonic) [9]. These terms will be 
written 


eo” = Lo Ci intitjhe , (i, j,k = 1, 2, 3), (17a) 


2. (17b) 


to 


eg’ = ye Dy 5,2 ;% E421 , (7, ds; ae 
where the indices run over the ranges given to yield no redundancy of symmetric terms 
(differing only in subscript order in the coefficients C;;, , D;;,,). Thus, the components 
of Ag’ and A ¢™ are homogeneous polynomials of degree two and three respectively 
in a; . But, from Eq. (10), any power x? (1 < gq < 3) can be expressed as a polynomial 
in s of lowest term s*. If the expressions derived in this way for x! are substituted in the 
expressions for VY ¢”’ and VY ¢“ corresponding to Eqs. (17a) and (17b) respectively, 
and powers of s above the third are dropped, one obtains 


if} =A+ Pir} + (E+ F)fs'}, (18) 


° . q 4) 
where the square matrices E and F are derived from the terms VY ¢'” and VY ¢™’ re- 
spectively of f in Eq. (16). Since VY ¢”’ is quadratic in z, , it can contribute no term 
to f which is linear in s; hence, E = [F,;] can be partitioned by columns thus 


E = (0, Ei2 , Eis). (19) 


Similarly, V ¢‘” can contribute no terms to f which are linear or quadratic in s; hence, 
F [F,,] can be partitioned 


F = (0, 0, F,3]. (20) 
The non-vanishing elements of E and F are given by 


Big = 8C,;078 + 2a; YS Ciiptp + D> CipghtyMe 5 (21a) 


E.=>p | 8c. + > C;,,,(a:8, + Bia,) + 5 D> Civalap8. + a,x.) | (21b) 


p.a.r 


Fis = AD i038 + 8a? YS Dy, ia) + 2a; SY Di ivetipte + D> DiperOtyQteXe (21¢) 


where the indices p, g, r run over the range 1, 2, 3 exclusive of 7 to yield no redundancy 
of terms differing only in order of subscripts in the coefficients C,;, , Dy jx: . Equation 
(18) contains the non-linear terms of V ¢ through order s* expressed as a linear function 


of s’ and s’°. 
Finally, by means of the linear relationship (13) between {s'} and {r}, Eq. (18) 


yields the quasi-linear approximation f~” to f, 
f°} =A+ (P4+Q)ir}, (22) 


where 


Q=(E+FS'r (23) 
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is a constant square matrix whose elements depend only on coefficients of g and the 
initial conditions of the motion (through the parameters a, 8, y, p, 7, p’). The parameters 
a, 3, y entering Q through I can be evaluated directly from standard formulas [8] as 


e= Vo 3= _Vo X (bo + A) - Xa, Y= a x 6B, (24) 


—_ kK 


Vo’ ' | Vo X (bo + A) | 


where v, (of magnitude v, ¥ 0) is the initial value of dr/dt, and bh) * —A is the initial 
value of b. The values of p, 7 and p’ in S™ are given by 


l/p = — | Vo X (bo + A) |, (25a) 
1/r = pto Vo X (bo + A)-do (25b) 
,’ = po [3Vo° (Do + A) — p$-d,], (25c) 


where the vector dy is defined by 
{do} = {(b)o} + Pf}. (26) 
The corresponding quasi-linear approximation to the differential equation (1) is 
ir} — (P+Q)irj = [bj +A. (27) 


The matrix Q vanishes and hence Eq. (27) reduces to the linearized equation (6) as 
p, tT — 0; thus the magnitude of the quasi-linear correction as well as the range of s over 
which it is valid is greater for p and 7 larger. 

The qudsi-linear equation (27) differs from the linearized equation (6) in an im- 
portant respect. Unlike P, the matrix P + Q is not symmetric in general, and thus 
Eq. (27) implies some unilateral coupling of its component equations. Since P + Q can 
have complex eigenvalues, reduction of the equation to normal coordinates is not 
generally possible. In view of this complication, analytic solution of the system (27) is 
best carried out by the Laplace transform method or other operational technique. 

The quasi-linear approximation (22) for f is equivalent to the expression of Eq. (18), 
which clearly contains all terms in f of order s* or lower. Since s = O(t) when v, ¥ 0, 
it follows that the Taylor expansion of the quasi-linear solution agrees with the true 
Taylor expansion of r through terms’ of order ¢’, as compared to the linearized solution, 
which agrees through terms of order ¢°. For v) non-vanishing, an integral curve of the 
quasi-linear equation has contact of fifth order with the actual integral curve at the 
initial point. The approximation is of osculating type in general, since it amounts to 
evaluating f on the osculating sphere of the integral curve at the initial point. 

The constraint that an integral curve lie in an equipotential surface requires a con- 
dition on b for a given f; the condition [10] is 


b+ f,. =(r)*/r, , (28) 


where b,, + f, is the component of b + f normal to the surface and r, is the radius of 
normal curvature of the surface in the direction of the curve. When the potential ¢ is 
spherically symmetric, the quasi-linear equation (27) becomes exact for integral curves 


*The partial sums through terms of order t§ of both Taylor expansions exist, because of the differ- 


entiability assumptions on \/¢ and b. 
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lying entirely in an equipotential. To show this result, select 2, in the direction of —V ¢ 
(of magnitude f,) at O, and let R, be the directed distance from O to the center of 
symmetry. From the symmetry of the field V ¢, one can write directly 


(r"} — folto'{r} = {b} + {0, 0, —fo} (29) 


for an integral curve on the equipotential ¢ = 0. For such an integral curve which is 
twisted at O, the partial sums through terms of order ¢° of the Taylor expansions of r 
are identical for the solution of Eq. (29) and the solution of the quasi-linear equation 
(27). This fact is sufficient to show that f,R,' is a triple eigenvalue of P + Q, and thus 


one has 
P + Q = foko'I (30) 


in this case, where I is a unit matrix. A corresponding theorem can be proved for the 
case where ¢ is cylindrically symmetric; the theorem for plane symmetry is trivial. In 
these special cases, the matrix P + Q is symmetric and all its eigenvalues are of one 
sign (even if g is harmonic). 

The preceding treatment is not valid for an integral curve which is plane (1/7 = 0) 
at the initial point O, because the inversion of Eq. (10) to yield Eq. (13) is not possible. 
In this case, if p be finite and non-vanishing at the initial point, an analogous treatment 
can be carried out by noting that S has a non-singular submatrix. In the corresponding 
development, of which the results will be stated, only terms in ¢ through the cubic 
> ean be included. Let T”* and A be the submatrices 


| 0 Qa, As As 
T= , A= , (31) 
0 2p B, B2 B, 
of S‘ and T respectively. If one writes 


R = (0, £,.]T"'A, (32) 


where the partitioned matrix [0, E.] is of 3 X 2 order, the corresponding quasi-linear 
approximation to Eq. (1) is 


fr} — (P+ R){r} {b} + A, (33) 


] 


which reduces to the linearized equation (6) for p — 0. For v) non-vanishing, an integral 
curve of Eq. (33) has contact of fourth order with the actual integral curve at the initial 
point, and the corresponding Taylor expansion of r agrees with the true Taylor ex- 
pansion through terms of order ¢*. 

4. Applications. The preceding results will be applied to the case where f represents 
a gravitational acceleration. The value of the coefficient matrix P of the linearized equa- 
tion which corresponds to a spheroidal earth rotating with angular velocity Q will be 
obtained from Eq. (8) on the basis of Clairaut’s first-order theory [11] of the variation 
of gravity. This matrix for the case of a non-rotating spherical earth is given in J. A 
solution of the quasi-linear equation (33) for plane motion corresponding to this value 
of P will be obtained. 

The initial point O will be taken on the surface of the earth (Fig. 1), with the 2;-axis 
directed along the outward normal to the surface, the z,-axis in the direction of a meri- 
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dian, and the z,-axis in the direction of a parallel of latitude (with right-hand sense). 
The potential g at points external to the earth can be written 
g=gat¢, (34) 
in which g,is due to Newtonian attraction, and 
@=5 (QF sin @) (35) 
is the potential due to centrifugal reaction at the point (R, 6), where R is the radial 
distance from the earth’s center and @ is the corresponding polar angle from the earth’s 
axis. The potential yg, can be represented as due to sources entirely within the earth 
and thus is harmonic at the earth’s surface; hence, the mass density o (which is fictitious 


in this case) is fixed by ¢, alone as 
= —-” (36) 


On the first-order theory [11], the potential (34) leads to a radius R of the earth for 
colatitude @, and a surface value g of the gravitational acceleration given by 
R = R.,(1 — A cos’ 8), (37a) 


q 


9 = Gel + B cos’ 6), (37b) 
where F,, and g,, are equatorial values of R and g respectively, and A and B are para- 
meters which depend on geometric constants of the earth, its mass, and its angular 
velocity. The lines of curvature of the (approximate) ellipsoid of revolution represented 
by Eq. (37a) are the meridians and parallels of latitude; the principal radii 7, and r, of 
normal curvature at O are equal respectively to the (negative) radius of curvature of 
the meridian at this point and the (negative) distance from this point to the polar axis 
along the surface normal [8]. Thus, one can write from Eq. (8), 


1+ Pi, 0 P%s 
P= —w,, 0 1+ P%, 0 : (38) 
51 0 ~~ P33 
where 
Wea = Jea/Rea 5 (39) 
and the assumption A, B < 1 yields 
” = 2A + (B — 3A) cos’ 8, (40a) 
P}. = (B — A) cos 8, (40b) 
2 2 ' 
433 = 2 + A+ = B + 2(B — 2A) cos 8, (40c) 
”, = Ph, = Asin 20. (40d) 
In determining Eq. (40c), use has been made of Clairaut’s theorem [11], 
5 ta 4 &, (41) 


Weg vo 
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so that the elements of the coefficient matrix P do not contain the earth’s angular velocity 
2 explicitly. 

To determine the coefficient matrix R of the quasi-linear equation for plane motion 
(1/r = 0) in a gravitational field, it is sufficiently accurate for the purpose at hand to 
take ¢ = g, R?/R, which corresponds to a non-rotating spherical earth of mean radius 
R, and corresponding mean surface gravity g, (as defined in J). The value of R is then 


R = w., [Ri], (42) 
where, for 7 = 1, 2, 3, 
Ri; = 6da,0;8;p/R., , (43a) 
Rj; = Gra.a38;p/R.q ; (43b) 
Rj; = 3d(1 — 3a3)8;p/R.. , (43¢) 
in which the factor 
N=14+2443B (44) 


arises from expressing g, and R, in terms of g,, and R,, respectively by means of Eqs. 
(37) and results in 7. Equations (42) and (43) generalize to an arbitrary plane of motion 
the corresponding result* in J. 

A solution of the quasi-linear equation (33) corresponding to P of Eq. (38) and R 
of Eq. (42) will be obtained for the special case of motion in a vertical plane at O. Take 


new axes z{ , x; at O (Fig. 1), where z{ is along the trace of this vertical plane in the z, , 


‘ 
X34 











Fic. 1. Coordinate frames on spheroidal earth. 


*The coefficient matrix Ap of Eq. (45) in J is in error as printed; the element in the first row and 
second column should have a negative sign prefixed. Otherwise, the results correspond with the identi- 
fication 2; , 22, %3 = 2X, Y, Z} a2 = Bz = Oin Eqs. (42) and (43) above. 
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z, plane and x; = 2; . Let 7 be the azimuth (positive sense from 2, to 22) of the 2{-axis 
relative to the 2z,-axis. The transform of Eq. (33) in the z{ , x3 frame is 


(xi"') 1+ Pn Dis (xi fi 

) + we, ( = ) , (45) 

3 Psi —2 — Psst {x3 b; — g, 

where 

Pi = Pi, cos’ 9 + Pi, sin’ yn — Ri, (46a) 
Ds3 = P33 + Ris (46b) 
Piz = Pi, cos 7 — Pls; (46¢e) 
Ds = Pj, cosy — Rj, , (46d) 


in which the R/, are defined by Eqs. (43) in terms of direction cosines @’, 8’ referred to 


the x/ , xj frame. Under the approximation | p;; | < 1, the solution of the quasi-linear 


system (45) corresponding to the initial velocity {v7 , vs} 1s 
r, = vw, sinw,t + « | bi(7) sin w,(t — 7) dr 
(47a) 
! “ 
+ 5 pisyvys(t) + | [bs(2) — g]s(t — 2) dre, 
o /0 } 
tu viws' sinh wt + w | [bi(r) — q| sinh w(t — 7) dr 
(47b) 
‘ at , 
+ 3 re | i8(t) + | bi(r)s(t — 7) dr, 
where 
| ) ”" | 
@, = w,1 + = Pui); @ = 2°‘a,.i 1 + 1 Ps f (48) 
and the function s(é) is defined by 
s(t) = w,;' sna,t — w; sinh wef. (49) 


linearized solution is obtained by letting p — 0 in the parameters 


The corresponding 
<< 1 imposed to obtain the solution (47) requires 


p,; . Note that the restriction | p,, 
p < R,, in Eqs. (43) for the elements of the matrix R. 

The author wishes to acknowledge helpful discussions of this problem with Dr. R. 
Isaacs, Dr. A. Latter, Dr. G. Peebles and Mr. H. Kahn of The RAND Corporation, 


and with Prof. T. Dantzig of the RAND consulting staff. 


Appendix 


In this appendix, Eq. (8) of the text for the matrix B will be derived. Take coordinate 
axes as specified by the text in connection with this equation. It is known that, in these 
coordinates, a plane parallel to the tangent plane of a surface and at the signed distance 
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p = x, from it cuts the surface in a curve (the Dupin indicatrix) defined within terms 
of second order by [8] 
2 2 
p = 1/27, + 22/2rz +--+. (ts =p), (1A) 
where r, and r, are the principal radii of normal curvature of the surface at O in the 
directions of x, and x, , respectively. From Eq. (2), this curve is likewise given by 


p= fi '(Byzy + 2B,.a2%2 + ByoX3) $+ -+e, (tz; = p), (2A) 
if f, = | V ¢| at O. Comparison of Eqs. (1A) and (2A) yields 
B,, = {fis B,, = 0, By, = : Sora’. (3A) 
2 2 
The coefficient B,, can be fixed by the condition that ¢g, in general, satisfies Poisson’s 
equation V° g = —e, where o is a mass density. This condition yields 
By = 5 (foKm + 07) (4A) 


where K,, is the mean curvature r,' + r,' of the -quipotential and o, is the mass 
density at O. 

To evaluate B,, , consider a plane normal to the z,-axis at the signed distance p = 2, 
from O. One can show that the equation of the curve of intersection of this plane with 
the equipotential ¢ = 0 is, within second-order terms [8], 

(1 + r2.:p ‘T2) = p, 2r, + v3 ret +s, (x, = P); (5A) 
where r,., is the derivative evaluated at O of the principal radius of normal curvature 
in the z-direction with respect to arc length on the line of curvature in the z,-direction. 
From Eqs. (2) and (3A), this curve is likewise defined by 


z,(1 — 2f5'Bisp) = p’/2r, + 22/2r2 + ---, (x, = p), (6A) 


which yields 
Row atte oe (7A) 
13 2 J0°2,1°2 
by comparison with Eq. (5A). Similarly, one has 
l -1 
Bas = 9 Sorat ) (8A) 


where r,. is the derivative evaluated at O of the principal radius of normal curvature 
in the z,-direction with respect to are length on the line of curvature in the z,-direction. 
These results vield Eq. (8) of the text. 
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SOME REMARKS ON A CLASS OF EIGENVALUE PROBLEMS WITH 
SPECIAL BOUNDARY CONDITIONS* 


BY 
G. W. MORGAN 
Brown University 
1. Introduction. ‘The study of a variety of physical problems frequently leads one to 
an eigenvalue problem. One very common eigenvalue problem, the Sturm-Liouville 
problem, consists of finding non-identically vanishing solutions of the differential equa- 
tion 
Lily(x)| + Awy = (py’)’ — qy + Awy = O, %SztE2, (1) 
subject to the homogeneous boundary conditions 
aoy’(to) — boy(to) = O, 
(2) 
ayy’(2,) + by(a,) = 0. 

Here p, g, w are functions of x which are positive in the domain and possess the 
required number of derivatives, \ is an unknown constant parameter, and do, by etc. 
are constants, such that ab, , a,b, > 0. 

In this paper we consider a problem which is identical with the Sturm-Liouville 
problem in all respects except for the appearance of the parameter \ in the boundary 
conditions. This introduces certain difficulties. Our primary object is to discuss a trans- 
formation which eliminates \ from the boundary conditions and which thus enables 
us to carry over several important results of the Sturm-Liouville theory. 

The method is then extended to higher order differential equations and finally 
physical examples are presented which provide an interpretation of the mathematical 
technique. 

2. Comparison with the Sturm-Liouville Problem. The Sturm-Liouville problem is 
known to have the following properties. 

Non-trivial solutions (eigenfunctions) exist only for certain positive values of A (the 
eigenvalues). The eigenvalues \, , A, «+ A; «++ , in order of magnitude, form a countably 
infinite set and the corresponding eigenfunctions y, , y. --* y, «+: constitute a complete, 
orthogonal system. Since the functions y, are determined up to an arbitrary constant 
multiplier only, they can be normalized. The orthogonality and normalization con- 
ditions are then expressed by 

" 0,i #J 
[ wy.y; dx = ) 
“Zo :. z = J 
By completeness we mean that any function f(z) which is piecewise continuous can be 
approximated as closely as desired in the mean square error sense by a sufficiently large 
number of terms of the finite series : ae c; y; , Where the coefficients c,; are given by 


(3) 


¢, = ‘% wfy; dz. (4) 


*Received April 30, 1952. 
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The problem we wish to study here consists of the differential equation (1) together 
with the following boundary conditions: 


y’ (xo) a K dy(Xo) = 0, 
y'(a,) — K,dAy(x,) = 0, 


where Ke, hs = 0. (5) 


We first note that the usual orthogonality relation (3) no longer holds. In fact it 
can be shown following standard procedures that 
| wy. y; dx = —K,p(x,)yi(a)yi(a) — Kop(ro)yi(to) yi(Xo), 1}. (6) 
For the Sturm-Liouville problem we know that a function f which satisfies certain 
conditions* can be expanded in a series 


f= Dew 
t=] 


where the expansion coefficients c; are given by (4). In our problem this is no longer 
true. 

Before discussing the transformation of the problem which will give, among other 
results, a physical interpretation of (6) and will also readily yield a formula for the 
expansion coefficients, we show how these coefficients can be derived formally by as- 
suming the validity of the expansion 


fxa=DdSKrvyi, mira. (7) 


Multiply by wy; and integrate over the domain, reversing the order of summation and 
integration: 


az, 


| why, dx = } is ¥: [ wy y; ax. 


/ £0 i=! “Ze 


Using (6) we obtain 


| why, dx = 7; | . wy; dx — p dg vy [Kip(ai)y(a.)y(2.) + Kop(xo)ys(Xo)y;(o) ] 


“2 i=] 


where >.’ indicates a summation in which the term 7 = 7 is omitted. To eliminate the 


coefficients y; for i ~ 7 we make use of the series (7) evaluated at the end points. 


f(%o) = viyil®o) + p by V:Y (20), 
s=1 


f(x) = Viy (21) + >’ veya). 


1=1 


_———————— 


*See for example Courant Hilbert, Methoden der Math. Physik, Vol. I, Ch. VI, Sect. 3.3. 
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We now have , 


ory 


| wfy, dx = 7; | wy? de — Keplae)ys(xo)(f(a0) — rsys(ae)] 


vr Zo 


— K,plx,)yi(a) (f(a) — viyi(x)). (8) 
If we require the y; to be normalized so that 
in wy; dx + Kop(xo)yi(xo) + Kip(a)yi(a) = 1, (9) 
then (8) gives 
+ * [ why; dx + Kop(xo)y;(to) f(t) + Kyp(a)yi(ai) f(a). (10) 


The above is essentially a generalization of the method employed by some authors in 
the special case when the y, are the trigonometric functions.* 

3. Transformation to Quasi-Sturm-Liouville Form. We now show how the problem 
can be transformed into, essentially, the standard Sturm-Liouville form. We shall here 
use heuristic arguments only, although it is possible to put the reasoning on a more 
rigorous basis by means of suitable limit processes. 

We define two functions 6(z, — 7) and 6(m — 2») in the following manner: 


d(x, — n) = O for n*#2, 


wail [ ” gil Mle > deed. 


| 
Similarly: 


i(n — X%) = 0 for ” Xo; 


and [ v(n) 5(n — ao) dn = v(2,). 
It is also convenient to define: 
af 0; 2s Ss gE<Xy 
h@) = | y(n) (x, — 1) dn = ; 
"oe ya); —E=Xy 
(11) 
se (0; Xo < é < Zi 
g(t) = | y(n) 6(n — 2) dn = ) 3 
™ y(t); § = Zo 


We now introduce a new dependent variable u by the relation: 


az pt 


u(x) = y(x) — AK, | | 


ZTo* Ze 


u(n) ae, ~ 9) dndg— Ko | u(x) &(y — 2) dn dt. (12) 


*Timoshenko; Vibration Problems in Engineering, p. 319. K. Ludwig; Einwirkungen von Bodene 
schwingungen auf Pfeiler, Z. a. M. M. 14, 361-363 (1934). 
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From (11) we see that 


| h(é) dé = 0, f= 25 fh. 


pe Bes (13) 


lA 


| g(é) dé = 0, Zo 


Hence (12) becomes u(x) = y(x),% <x ay. 
We now find the differential equation and boundary conditions satisfied by the 


function u(2). 
From equations (12) and (11) we have: 


u(x) = y'(2) \K, | y(n) 6(a, — n) dn + AKo | y(n) 6(n — 2) dn, (14) 


“2 


or u'(x) = y’(x) — AK, A(x) + AKog(z), 


and therefore 
u(x) = y'(2), <2 <i mh, 


u'(a,) = y’(x,) — AKyy(x), (15) 
u'(Xo) = y’(Xo) + AKoy(Ao). 
Comparing with (5), the boundary conditions on y, we see that our boundary con- 
ditions on u are 
u'(z,) = u'(z,) = 0. (16) 
From (14) and (13) we obtain: 


A''(x) = y’"(x) — AK, G(x, — x)y(x) — AKo G(x — 2) y(z) 


— 2)u(z). 


~ 


= y’"(x) — AK, 6a, — x)u(xz) — AKo 
Substitution into the differential equation (1) yields 
(pu’)’ — qu+ Alwu + Kop h(x — x)u(x) + Kip 6a, — x)u(z) 
— Kop’g(x) + K,p’h(x)] = 0. 


Now g and h vanish everywhere except at the end points, and at the end points they 
are irrelevant in view of the presence of the terms containing 6(z — zp») and 6(2, — 2), 
respectively. Hence the differential equation for u can be written: 


(pu’)’ — qu+ Apu = 0 (17) 


where p = w + Kop(2) 6(x — Xo) + Kip(x,) 6(@, — 2). 

The problem defined by (16) and (17) is in the standard Sturm-Liouville form except 
for the singularities of the function p at the ends. We shall formally regard it as a Sturm- 
Liouville problem and assume without attempting a rigorous justification that many 
results of the Sturm-Liouville theory, such as those concerning the nature of the eigen- 
values, the completeness of the eigenfunctions, etc. are still applicable. 
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The orthogonality (3) yields 
o - 0,0 #J 
| puu,; dx = | py.iy; dx = * . 
es by 1 = J 


Substituting for p(x), the integral becomes 


| wy .y; dx + Kop(2) | yy; 1x — Xo) dx + Ky,p(z,) | yy; (a, — 2) dz, 


and this gives 
si 0,0 #) 


| wy.y; dx + Kop(to)yi(Xo)y;(%o) + Ky p(a)yfa)y(a) = (18) 


t= 
This is immediately seen to be identical with (6) and the previously introduced 


normalization (9). 
The coefficients y; in the approximation 


n 


f sinc > yu; = dvs 


i=l 


are immediately found by using the appropriate functions in formula (4); 


az, ara 


n= | pfu, dx = | pfYyn dx 


“Zo To 


| why, dx + Koplae)f(to)ya(te) + Kipla.)f(x)ya(as), 


“z 


II 


and this is identical with (10). 
Completely analogous results can of course be derived if \ only appears in one of the 


boundary conditions, the other condition being as in (2). In that case only one of the two 
double integrals in (12) will have to be used in the definition of the new variable u(z). 
4. Extension to Higher Order Self-Adjoint Differential Equations. An analogous 
method applies to the equation 
(poy’’)’’ — (py’)’ + pu — Awu = 0, % StS. 
This is to be solved subject to four homogeneous boundary conditions, one or two of 


(91) 


which involve the parameter \ in the following manner: 
y’’"(xo) — Kody(%o) = 9, 
y’’"(x,) + Kiry(x) = 9, (20) 
K, , K, > 0. 


As an illustration consider the case in which, in addition to (20), the conditions 


y’’ (ao) = y’’(x,) = 0 are to be satisfied. 
Define the new variable u(x) by 


ar ar as 


u(x) = y(x) + AK, [ 


“ro (to ©20o Y2 


+—Kof ff / y(n) 8(n — 2) dn de ds dr. (21) 


“Zz 


F 
y(n) 5(x, — ») dn dé ds dr 
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Proceeding as before we obtain 


u(x) = y(2), 


u(x) = y’(2), (22) 
u''(x2) = y’'(2), oo = F = Se, 
Also 
ul (a2) = y!'""(x) + AK, A(x) — AKog(x) 
where g and h are defined as before by (11). Hence: 
u’’"(z) = y'’"(2), x2 <2, 
ula) = y’""(a1) + AK, y(a)), (23) 


u(t) = y'"" (to) — AKoy(Xo)- 
Comparing (22) and (23) with (20) we see that the boundary conditions on u are: 
u(t) = ul'"(a) = u''(z,) = u’""(z,) = O. (24) 
Finally we have: 
u'Y = y'* + XK, 8&2, — z)y(x) + AKo S(x — xo)y(z). 
Substituting in the differential equation (19), we obtain 
(pou’’)’’ — (p,u’)’ + pu — Apu = 0, (25) 


where p = w + Kopo(Xo) 6(@ — Xo) + Kip2(x,) 6(a, — 2). 

Hence we have again succeeded in converting the original problem into a standard 
form, the differential equation having been changed only by the replacement of the 
original weight function w by a new weight function p identical with w everywhere 
except for singularities at the boundaries. We are again in a position to apply the results 
of the standard theory. 

The orthogonality relations and the formulas for the coefficients 7; in the approxima- 
tion 

n ‘ 
f~ Dru = Dry: 
s=1 s=1 
are found precisely as before. 

This technique can be extended in a self-evident manner to self-adjoint differential 
equations of order 2m, m = 1, 2, --- if one or two of the boundary conditions are of the 
form 

y?""(z,) + KaAy(2,) = 0, i=0,1 


and the K; have the proper sign to make p(z) positive. 

5. Physical Interpretation. We now consider three illustrative examples of physical 
systems, the mathematical analysis of which leads to the type of eigenvalue problem 
under discussion. 

(a) Longitudinal vibrations of a bar with a mass at the end. 
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The differential equation is 
av 4] ( 2) 
w—3 = — | AE — 26 
ot Ox Ox (26) 
where v is the longitudinal displacement of a cross section of the bar, w is the mass per 
unit length, A is the cross section, and E is Young’s modulus. 
If the bar is built in at z = 0, the boundary condition there is v(0,t) = 0. If a mass 
M is attached to the end x = I, the force in the bar at that end must be equal to the 
rate of change of momentum of the mass M. 
Hence the boundary condition at « = 1 is: 
, Ov av 
AE — = —M = at z=l, t> @. (27) 
Ox dl 
To solve this problem we look for normal modes in the form 
v(x, t) = sin wt} 
-y (x) 
cos wt) 
and hence obtain the following differential equation and boundary conditions for y: 


(AEy’)’ + AXwy = 0, 
y= 0 at x = 0, 
AEy’ — My = 0 at z=l 


where = w, | (28) 


This is precisely of the form (1) and (5) with a standard boundary condition at z = 0, 
Transformation of this problem’‘to the form (16) and (17) gives: 


(AEu’)’ + (w + M O(l — z))u = 0, 
u=0 at z=0, 


u’ = 0 at z= l, (29) 


The physical meaning of the formulation (29) is that the mass M is treated as a part of 
the elastic system by making the mass per unit length infinite at z = / in such a manner 
that the total mass added is M. The new boundary condition at z = / is now applied 
at the outer end of the mass M and is the condition for a free end. 

The new orthogonality condition, given by (18) is: 


al 


[ (w+ M ol — x)lyy; dx = | wy .y; dx + My(Dy() = 0, i). (30) 


“0 


It clearly shows the physical meaning of the transformation. The infinite mass density 
at « = | changes the weight function to account for the extra mass /. 
(b) Lateral vibrations of a beam with end load.* 





*F. Gassmann, Uber Querschwingungen eines Stabes mit Einzelmasse, Ingenieur—Archiv 2, 222-227, 


1931. 
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The mathematical formulation of this problem is: 
dv a ( , a) ; 
wW—> ;\ HI —3}) = 0 (31) 
ot + Ox” Ox” 


where v is the lateral displacement of a cross-section and /7/ is the flexural rigidity. 
If the beam carries a mass M at x = I, the boundary conditions there are: 


xn OF oy 
EI — = M —3 
Ox ol 
and 
Ov 
——_ am (. (32) 
Ox 


These conditions assume that the center of gravity of the mass M is very close to 
x = /so that terms arising from this distance can be neglected. 
The search for normal modes leads to: 


(Ely’’)’”’ — rXwy = 0, (33) 
and the boundary conditions at x = 1 become: 


Ely’” + \My 


I 


(34) 


vr 


yy’ = 0. 


This is of the form discussed in section 4. If we carry out the transformation we again 
see that its physical interpretation is the inclusion of the mass M into the elastic system 
with the effect of obtaining a new weight function w + M é(l — 2). 

(c) Heat Conduction in a rod with finite energy reservoir at end. 

Consider one end of the rod to be in contact with a body in which heat diffusion is 
so rapid that it is permissible to consider the body to be at a uniform temperature equal 
to that at the point of contact. We also assume that the only energy flux into, or out of, 
the body occurs at the contact with the rod. 

We have 

moi = az (KA 5a) (35) 
where v is the temperature, AK the termal conductivity, c the specific heat per unit mass, 
m the mass of the rod per unit length. If, at x = /, the rod is in contact with a body of 
mass M and specific heat y, then the boundary condition there is: 


Ov Ov 
—Ks = 7M —, z= l. 36 
KA rs yM al x (36) 
We try solutions: 
v =e ™y(z). 


The equation for y is: 
(KAy’)’ + XA(cm)y = 0. (37) 
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The boundary condition at z = 1 is: 


KAy’ — d\yMy = 0. (38) 
Comparing with (1) and (5) we see that 
KA =p, cm = w, and yM/KA = K,. (39) 


Our transformation yields: 


(KAu’)’ + Mem + yM Sl 


x))u = 0, 


and 
u’ = 0 at z= l, (40) 
The physical interpretation is the same as in the previous examples. The body at 
the end is thought of as a part of the system in which conduction takes place. This is 
done by making the heat capacity per unit length infinite at x = 7 in such a manner 


that a total heat capacity y.V/ is added to the system at that end. 
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THE CHEBYSHEV APPROXIMATION METHOD* 


BY 
PRESTON R. CLEMENT** 


Electrical Engineering Department, Princeton University 
g g I 


A treatment of the general Chebyshev approximation method as it interests physicists 
and engineers is given, with a detailed discussion of the properties of Chebyshev poly- 
nomials. Brief applications to electric circuit theory are presented.*** 

Introduction. Of the various means of approximating a given function, the Cheby- 
shev method is, for physicists and engineers, one of the most interesting and important. 
Originally introduced by P. L. Chebyshev in a mechanical linkage problem, (26) this 
procedure came into particular importance in electrical engineering with the publication 
of a new method of filter design by W. Cauer.(5) 

It is the aim of this paper to review the various articles written on this subject, with 
the purpose of presenting the material in a manner sufficiently general to permit physicists 
and engineers to appreciate fully this approximation problem. Particular emphasis is 
given to Chebyshev polynomials, with brief applications to electric circuit theory. 

1. The Problem of Approximation. Chebyshev Approximation. Consider a function 
f(x) defined in an interval a S x S b. Let g,(zx; a;) be a function defined in the same 
interval with n adjustable parameters a, , --- , a, . Our problem is the following: What 
should be the values of the n parameters in order that g,(x; a;) will approximate f(x) as 
well as possible in the interval (a, b)? 

The answer to this question depends, of course, on what is meant by the expression 
“as well as possible.”” The applied mathematician is interested in obtaining a function 
J,(2; a;) which will be easier to use in his calculations than the original function f(z), 
yet which will produce a result sufficiently close to the exact answer. A frequently used 
approximation is that of least squares; that is, the a’s are adjusted so that 


| (f(x) — g(x; a,)]? dx 
is a minimum. For example, if the interval (a, 6) is (—z, x), if f(x) is odd about the 
origin, and if g,(x; a;) = a, sin x + a, sin 2x + +--+ + a, sin nz, then the least squares 
approximation will produce the Fourier coefficients for the a’s. 

Let us now consider an engineering problem, the design of filter networks on an 
amplitude spectrum basis. Suppose that it is desired to have a transmission characteristic 
f(x) which has a constant non-zero value in the pass bands and a zero value everywhere 
else. It is clear that with a physically realizable circuit such a characteristic can only be 
approximated. Thus in this case, we can denote by g,(2z; a;) the amplitude spectrum of 
the actual filter circuit which should approximate our ideal characteristic f(x) “‘as well 
as possible.” x will denote the frequency, and the a, will be numbers determined by the 
adjustable circuit parameters of the system. 

For the purpose of filtering, it is clear that the least squares approximation may not 
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be a very satisfactory one, since one may desire the guarantee that in the stopbands the 
attenuation for every frequency is above a certain definite value. With the least squares 
approximation, the attenuation in the mean of a frequency mixture in the stopbands is 
indeed very high, but for one or more frequencies it may be quite small. 

Thus to insure that such large departures from the ideal characteristic do not occur 
for any frequency, we would demand that the maximum of 


| f(z) — g,(x; a,) 


be a minimum in the prescribed interval (a, b). When such is the case, g,,(x; a;) is said to 
approximate f(x) in the Chebyshev sense in the interval (a, b). Hereafter, we will mean 
by “approximate,’”’ approximate in the Chebyshev sense. 


If we let h,(x; a;) = f(x) — g,(x; a;), our general problem becomes the following: 
Given an arbitrary function h,(x; a;) with n parameters a, , «++ , a, , what should be the 


values of these parameters to make the maximum difference from zero of h,(x; a,) a minimum 
fora = x S b;i.e., so that h,(x; a;) approximates zero in (a, b)? 

2. The Function },(x;a;). In all that follows in this section, it will be assumed that 
h,(z; a;) is an arbitrary function defined in the interval a S xz S b, subject to the re- 
striction that the function and its derivatives with respect to x, a, , --: , a, are finite and 
continuous. In order to determine the values of the parameters for the maximum devia- 
tion of h,,(z; a;) from zero to be a minimum, the actual value of this deviation, and the 
points at which this deviation is assumed, two sets of equations will be considered. 

The first set arises simply by noting that if L, is the smallest possible maximum 
deviation of h,,(x; a;) from zero in (a, b), then when the parameters are properly adjusted, 


the equation 
h(x;a;) — L, =0 


will give the set of points x, , 22, , 2, in (a, b) at which the maximum is assumed. 
If the value L,, is assumed on the boundary of the interval, then two of the points of the 
set z, , ::: , 2, are a and b. The rest of the points are on the interior, where the first 
derivative must vanish.Hence if differentiatiun with respect to z is denoted by a prime, 


the first set of equations to be satisfied is 
W(z:a,) —-L? =0 ( 
(1) 


(2 — a)(x — D)hi(x;a;) = 0.4 


These equations will have « common solutions which are real and distinct; x, , 72, °°: , 
Zp. 
The second set of equations comes from a theorem stating a necessary condition on 
h,(2z; a;) in order for the function to satisfy the Chebyshev condition. 


Theorem 1. Given a function h,(z; a;) continuous along with its derivatives in the 


] 


interval a < x S b. Let 2, , --- , x, denote the set of points in (a, b) at which | h,(z; a,) | 
assumes a maximum. Let the greatest of these maxima be denoted by M. Then if M = 
M(a, , a2, *** , @,) has a minimum value for the parameters specified (i.e., M = L,), 


the following set of equations is insoluble for every set of non-zero real numbers S, , S2 , 
--+ , 8, such that the sign of S, is the same as that of h,(a,; a;),k = 1,2, --- ,u. 
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Ih, (2, ; a) Oh, (ey 5 a, Oh, (21 5 a, 

Oh, a a ‘‘. 4. Halas a;) 1 wee 4 Arnltr 54s) = §, 
0a, 0a, 0a, 

Dale 54s)» 4 Dhalte sa), 4... 4 Walt O), _ (2) 
0a, 0a, 0a, 

Dh (Xy ; a) - Oh, (Zp 34) ae Serre Ah(Xy 5 as) = S, 
Oa, 0a» 0a,, 





A theorem similar to this was proved by Chebyshev (27). The extension to the present 
form was given by Kircherberger (20). 


Example 1. h,(x; a;) = x? + ax + a,,(-1 S$ 2 S 1) 
For this function, equations (1) become: 
(x? + a,x + a.) — L; = 0, (3) 
(x — 1)(a + 1)(2x + a,) = 0. (4) 
Substituting two solutions of (4), z = + land z = —1, into (3) and subtracting the two 
resulting equations: 
a, + a,a, = 0. (5) 


Taking a, equal to —1, which certainly satisfies (5), the desired polynomial appears to be 
a + 2r — 1 (6) 


with a maximum deviation in the interval (—1, 1) of 2. 
Let us use the plus sign in (6). Then the polynomial has a maximum at + 1 of + 2 
and a minimum at —1 of —2. Equations (2) then become 


—-A, +r. = -S, 


Ai + Az S, 
where S, , S, are positive. Thus 


A = (S: + S,)/2, 


A, = (S, 7: S,) ‘2, 


violating the necessary condition for Chebyshev approximation given by theorem 1 
that equations (2) be insoluble. The same would apply if the negative sign had been 
chosen in (6). Thus (6) is not the desired polynomial. 

To solve our problem we note that another solution of (5) is possible, namely, a, = 0. 
Substitution of this value for a, gives for a, the value of —4}. Hence the polynomial 


becomes 


4 (7) 
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with a maximum deviation of } at x = + 1, 0. For this function equations (2) become 


—A, + rz — Si, 


Ao = —S, (8) 


? 


A +A. =8;, 
where S, , S, , S; are greater than zero. It is clear that (8) is insoluble in the \’s. Thus 
the necessary condition that the polynomial approximate zero in the Chebyshev sense 
is satisfied. That this is the proper polynomial of second degree follows from the fact that 
the maximum deviations are equal and three in number. This property of polynomials 
will be shown in section ITI. 


Example 2. h,(x;a,) = x + a, sin (72/2), (—-r S22 Sr). 


Again the problem is to find the value of a, such that the maximum difference of 
this function from zero in the interval (—7, 7) is a minimum. 
Equations (1) becomes 


x + 2a,x sin (77/2) + aj sin’ (72/2) — L; = 0, (9) 


(x — r)(x + w)[1 + (7/2)a, cos (7x/2)] = 0. (10) 


By a process similar to that used in example 1, it can be shown that the function z + a, 
sin (7x/2) which satisfies the Chebyshev condition in the interval (—7, 7) is 


h(x; a,) = x + 0.39 sin (72/2), 


with its maximum deviation equal to 2.75. Its graph is shown in Fig. 1. 
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Fic. 1. The function z + a; sin (7z/2) which satisfies the Chebyshev condition in the interval (—2,2). 
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Example 2 was presented primarily as an illustration of a function which satisfies 
the Chebyshev condition yet whose relative maxima are not all equal. For many general 
functions, however, the Chebyshev condition is satisfied only when all the relative 
maxima are equal. This statement is made more precise in the following theorem. 


Theorem 2.* Given a function h,(x; a;) defined and continuous along with all its 
derivatives in the interval (a, 6). Let « be the number of ordinates of h,(x; a;) where the 
magnitude of the function is a relative maximum, including the end points a and 6 and 
let « = n + 1, where n is the number of parameters. It is further assumed that either 
the first derivatives of h,(2z; a;) with respect to the a; do not vanish at z, , --* , 2, or 
else if the first derivative vanishes, the first non-vanishing derivative is of odd order. 
Then if h,(x; a,) satisfies the Chebyshev condition in (a, b) the y relative maxima are all 
equal. 

This so-called equal maximum property of Chebyshev functions will also occur in 
certain special cases where » # n + 1. For instance, in example 2 if the antisymmetrical 
function had been taken as x + a, sin (52/2) in the interval (—7z, 7), then all four 
maximum deviations would be equal under the Chebyshev condition, even though there 
is only one parameter to adjust. 

We will now leave the general Chebyshev function and consider the important 
special case of the Chebyshev polynomials. 

3. Chebyshev Polynomials. Given a function f(x) which is finite and continuous 
along with all its derivatives in the interval (a, b). Let g,(z; a;) = a,x" * + az"? + 

+ a, be the approximating function. Hence the function which we desire to approxi- 


mate zero in the Chebyshev sense is 

h(x; a,) = f(x) — fav" + az"? + --- +,]. (11) 
Let 2; , 22, °*: , x, be the set of points in the interval (a, b) at which the function given 
by (11) assumes its maximum deviation L, . Under these conditions, equations (2) 
become 


itt + Aa? + ees AR +A = S,| 
Asa + Aor? ‘ + ++ + Anite FA, 


II 
o, 
nN 


(12) 





Az, + Aor, . + a + Ren125 + An _ S, 


Now let P(x) denote \,2""* + \.@"-? + --- +A,, an entire rational function of (n — 1)th 
degree, which by system (12) is determined only in so far as that at a prescribed point it 
has a prescribed sign. To settle the question of whether equations (12) can be satisfied 
by a suitable choice of coefficients, it is necessary to consider the number of changes of 
sign of S,, S.,---,S,. If there is a change of sign in going from S, to S,.,, , then there 
is a root of P(x) = 0 between x, and z,,, . Since we have (n — 1) roots at our disposal, 
and since the number of changes in sign can never exceed the number of roots, then it 
follows that if the number of sign changes in the S’s is less than or equal to (n — 1), 


*The statement of this theorem and its proof follows along the lines of one given by LeCorbeiller (22), 
except that it was found necessary to add a hypothesis to take care of special functions. 
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system (12) will be soluble for the \’s. Hence the following proposition has been demon- 


strated; 
A necessary condition for the function h, (a; a;) as given by (11) to approximate zero 


in the Chebyshev sense is that among the set of values h,(z, ; a;), Aa(%2 3 a:), +++, An 


(x, ; a;) there are at least n changes in sign. 

From this it follows that if h,(z; a,) as given by (11) is to approximate zero, then 

equation set (1) will have at least (xn + 1) common solutions in (a, 6), and further the 

maximum deviation L,, will be alternately a maximum and a minimum of the function. 
To show that the condition of the non-solubility of equations (2) is not only necessary 


but also sufficient for the function given by (11), it will be shown that the condition 


leads toa unique function h,(x; a;). 


To prove this, assume that there exists a function /,(2; a,) = f(x) — (aa" + a 
r"* +--+ + a,) which also satisfies the non-solubility condition, in other words, whose 
maximum deviation is less or equal L,. Thus it follows that (h, — h,) is a polynomial 
of degree at most (n — 1). By the previous discussion there are at least (n + 1) points 
at which | h, | = L, and at least (n + 1) points at which | /, | < L, . Thus (h, — h,) has 
at least n zero points in (a, b). This means that necessarily h, = h,, . 


Let us now take for {(x) the simplest possible form: 


f(x) = 2" (13) 
Thus the function which we wish to approximate zero may be written as 
h(a: b,) = 2" + ba" + baw" + ee +O, (14) 
where the b’s of equation (14) are simply the negatives of the corresponding a’s of equa- 
tion (11) 
For convenience, let us repeat equations (1). 
he (x: b,) - | br = (), (la) 
(x — a)(x — b)hi(x; b;) = O. (1b) 


Since equation (1b) is of degree (n + 1), it follows from the previous discussion that all 
its roots must be distinct, and furthermore they must all also be roots of equation (1a). 


Now since the roots of A’(x: b,) = 0 are roots of (la), it follows that there are (n — 1) 


double roots of (la). Hence the system of equations 


2 


hi(z; b;) -— L, = 0 


[hi(a; b,)}’ = (0) 


will have (2n — 2) common solutions in the form of pairs of double roots. The other 
two common solutions to (la) and (1b) must occur at x = a and x = b. Thus the equations 


h(x; b;) — L, = 0 


(x — a)(x — D)[AK(2; b,)} = 0 


will have 2n common solutions. When it is recognized that the coefficient of 2°" in the 


first of these equations is unity, while in the second it is n’, it is seen that the differential 
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equation to be solved is 


oe , 
i-2= 2 (ate) (15) 
dx 


n 
There are four solutions to equation (15) which result from the four possible combi- 
nations of signs appearing after the square root of both sides has been taken. These 


solutions are 


sin f fsin\ —2r +at+b | 
h =. n are ( + >? (16) 
cos | leos/ a—b J 
where C is the constant of integration. Since our problem is to determine the polynomial 
of nth degree which approximates zero best, we are restricted to the solution 


—27r I - 
h, = L. cos E are cos ( Seb & ) + c}. (17) 


a— b 
The constant of integration is determined by the requirement that at the boundaries 
of the interval (a, 6), the polynomial must assume its maximum deviation. This re- 
quirement is satisfied for C = 0, though, of course, placing C equal to an integral multiple 
of z would not change the function A, . 


17) in a form which shows its polynomial character more readily, let 


—2r + a +) 


In order to get 


2 = arc cos ( 


a—b 
Then 
h, = L, cosnz = L,/2|(cosz + isin z)" + (cosz — isin z)"] 
= 5 (-2e tat b+ Via" — 2a + De + dab)” 
“aa — 0 


+(-2e7+a+b— V42" — Aa+ bx + 4adb)"]. (18) 


For the coefficient of x" to be unity, as given in (14), 


(a — b)" 
LL. = a a (19) 


Tabulation of Chebyshev polynomials in the interval (a, b). 
From equation (18) the following set of polynomials are obtained; 


Sa ew a us bh 
h, = om 3(a + b) i i a’ + Gab + h? 
y S 
» __ 9a + db) 3a’ b + 3b? 7 
fi ab ng EG teat £3 r— 1S? @ + Mad + 8) 


Chebyshev polynomials in the interval (—1, 1). 
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It is usually convenient to deal with the Chebyshev polynomials in the interval 
(—1, 1) instead of the interval (a, b). From (17) and (19), these functions may compactly 


be written as 


] 
— COs (m arc COS 2). 
9-1 


For our purposes the polynomials of unit maximum amplitude in the interval (—1, 1) 
will be more useful: i.e., the polynomials given by 
7’,.(z) = cos (n are COs 2). (20) 


Hereafter, references to Chebyshev polynomials will means those functions given by 
(20). The first ten such functions are tabulated in Table I for convenience. Fig. 2 shows 


a few of these polynomials. 


TABLE I. TABULATION OF THE CHEBYSHEV POLYNOMIALS 


T,(z) = cos (n are cos 2) 


T(x) = 7 


T(x) = 82° — 82° + 1 

T(x) = 162° — 202° + 52 

T(x) = 322° — 482° + 182° — 1 

T(x) = 642° — 1122” + 562° — 7x 

T(x) = 1282° — 2562° + 1602° — 322° + 1 

T(x) = 256a° — 5762" + 4322” — 120z° + 9x 

T(x) = 512z2'° — 12802° + 1120z2° — 400x* + 502” — | 


Another convenient form for 7',(x) is, from (18), 


T(x) = lig + V2? — 1)" + (2 — V2’? — 1)’). (21) 
-_ 


Also, the expansion of (20) gives: 
= n—A n — 
T(x) = (—( ) eee <eoeeeentese (n=1 (22) 


where [n/2] denotes the largest integer contained in the set 0, --- , n/2. 
A still further form for 7,,(z), which may be checked by using the original differential 


equation, is 


! aan eo = coe 7 Qn-1/5 (oO. 
T(z) = 2"(—1)" nyt V1i—2z = a-sy"* 29 (23) 
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Further properties of the Chebyshev polynomials. 1. Recurrence relation. Letting again 
z = are cos 2, T,, = cos nz. Since cos (n + 1)z + cos (n — 1)z = 2 cos nz cos 2, then 
T..; + T,-, = 27,7, . Thus taking into account from (20) the values of T,, for n = 0, 1: 


Tia, = 22T, — T,-1 (24) 


where 





“ 
a 











Fic. 2. The Chebyshev polynomials 7,,(x) = cos (n arc cos x). 


Ie fe 


2.Orthogonal property. The 7’s are othogonal to each other in the interval (—1, 1) 
with respect to the function 1/-V1 — 2", as is easily seen from the fact that 
0 (m # n) 


Jw 


2 


II 
IV 


| T.(2)T.(2) : p-— 5 = [ cos n@ cos m6 dé = (m 1) (25) 


nn Ving ; 





\xr(m =n = 0) 


3. Generating functions. One generating function which may be used to define the 


T',(x) is: 


1 — zt ‘it ‘ 
i-m+e?" p> Tat.  (t< 1) 


A second such function is 





I — ¢ 
h ——_—_—— = T(x) —. <1 
VFoaaei tO <D 


4. Maximization of the largest root. This property will be stated in the form of a 
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2 direction, as shown in Fig. 4. Thus we may say that the curve is given in space by z = 
cos n@ and r = 1. Therefore, the abscissa of the orthogonal projection of this curve on 
the plane formed by the z-axis and the @ = 0 line is then x = cos 6. Hence the projected 

















Fic. 4. The curve cos 76 on a cylinder of unit radius and its projection on a plane to form a Lissajous 


nattern. 


figure gives the relation between x = cos @ and z = cos né@ and is thus represented by 


= COS (7 arc COS 2). 


» 
~ 


If time is introduced as the variable instead of 6, the projection of the curve z = cos 


n@ can be written in the parametric form 


z= cos nt, 


x = cost, 


and thus represent Lissajous figures for an integral ratio of the two composing frequencies. 

If the curve had been taken as z = sin n@, then the projections would have been 
proportional to the U,,(x) functions instead of to the 7’,(z). 

9. T,(x) and U,(z) in relation to Fourier series. A similar geometric interpretation 
can be given (30) to show the relation of the 7, and U,, functions to Fourier series. The 
connection can be seen, however, quite easily in the following way. Let Y(6) be a given 
function in an interval, say (a, b). Developing this function in a Fourier series, we have 


Y(6é) = >. (a, cosn@ + b, sin n6). (as 6c Db) 
If x is set equal to cos 6: 


Y(cos' 2) = >> [a,7,(z) + b,U,(2)]. (a S$ cos’ zx S b) 
10. Expansion of an arbitrary function in Chebyshev polynomials. By taking into 
account the orthogonality property, it is an easy matter to show that if the series 


2c) + >. c,T,( x) 


n=1 
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converges uniformly to the sum f(z) in the interval — 1 S x < 1, then 


1 
«=? s@)TA) —4—. (30) 
wJ-1 V1 — 2x 
As an example of the use of a function expanded in Chebyshev polynomials, consider 
the curve of plate current versus grid voltage for a triode. (31) Linear theory of triode 
amplification is based on the assumption that the operating region is so small that the 
characteristic may be replaced by a straight line throughout the region with negligible 
error. This involves considering only the linear term in the power series expansion about 
the operating point. To calculate the harmonic distortion resulting from a finite grid 
swing, the higher derivative terms in the power series become important. From experi- 
mental measurement the coefficients of the Taylor series can only be derived from the 
limit towards which the measurements converge for an infinitely small grid swing. 
Consider a sine wave applied to the grid, so that the grid voltage e,(¢) is given by 


e,(t) = €0 + a cos at. (31) 
The application of e,(¢) causes a plate current to flow of the form 
i,(t) = In + 1, cos wt + I, cos Qwt + +++. 
Now if we write 
x= coswt (32) 


then 


7,(t) _— Io + I,T,(2) + I,T (2) + e**. (33) 


In other words, the amplitudes of the higher harmonics which are measured ex- 
perimentally, are, for a given operating point and grid swing, the coefficients of the 
triode characteristic expanded in a series of Chebyshev polynomials. Thus the expression 
of the triode characteristic in this form is easy to obtain experimentally, and when once 
written, gives the information regarding the amount of harmonics under particular 
operating conditions. 

Application of Chebyshev Polynomials to Antennas. It is frequently the problem in 
antenna design that the beam should be as narrow as possible, the power gain should 
be a maximum, and the sidelobes should be relatively small. These requirements are 
normally not easy or even possible to satisfy simultaneously, so the problem arises as to 
what to term the optimum field pattern. For this discussion we shall term that field 
pattern optimum which has a minimum beam width for a specified sidelobe level. 

To determine a method for obtaining the optimum pattern, consider an array of n 
isotropic point sources of uniform spacing and all in the same phase, as shown in fig. 5. 
For even n, the array is given by 5a; for odd n, by 5b. The amplitudes of the individual 
sources are denoted by the A’s as shown. Hence the total field for 5a in a direction @ is 
given, for large distances, by the expression (21) 


E..en = 2A, cos 6/2 + 2A, cos 36/2 + +++ + 2Ay,/2 cos (n — 1)6/2, (34) 
and for 5b, 
E.aa = 2Ay + 2A, cos ¢ + 2A, cos 26 + eee + 2A (ani) /2 cos (n a 1)¢/2, (35) 


where @ = 2rd/) sin @ and d is the distance between point sources. It is clear that (34) 
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Fia. 5. Line al broadside array of n isotro vie sources with uniform 8 yacin ° 
i I 


a. niseven. b. n is odd. 


and (35) may be written in the form: 


E..en = 2 >, A, cos (2k — 1)6/2 = 2 D> AT n-,(2), (36) 
k=1 k=1 

E.u=2 > Ajookp=2 > A,T.A(2), (37) 
k=0 k=0 


where z = cos ¢/2. 

It is now the problem to determine the value of the A’s in (36) and (37) so that the 
optimum pattern is obtained. One way is to make all of the A’s equal. This, indeed, gives 
a Maximum gain and a narrow beam width, but the level of the side lobes is very high. 
Another way is to make the A’s proportional to the binomial coefficients, in which case 
the side lobes disappear altogether if the spacing between elements is less than half a 
wave length, but unfortunately increased beam width and loss of gain result. 

Now let us consider Dolph’s method (14) of using the Chebyshev polynomials to obtain 
an optimum pattern. Suppose we have an array of n sources as in Fig. 5, where n may be 
odd or even. Let the specified ratio of the main lobe maximum to the minor lobe level 
be R. The problem is to find the relative source strengths for the pattern to be optimum 
in our prescribed sense. 

The forms (36) and (37) show that the field due to the n given sources is described 
by a polynomial of (n — 1)th degree, which may be denoted by P,-_,(z). The solution of 
the equation P,_,(z) = R will give zg , the point at which the pattern assumes its 
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maximum value. Letting y = 2/z, , then (1/R) P,-,(y) will have a value of unity at 
y = 1, and we have that y (rather than z) is equal to cos ¢/2. It is necessary that P,_,(y) 
have all its roots in the interval (— 1, 1), and furthermore be an even function about the 
origin. Thus P,_,(y) satisfies the requirements of Theorem 3. Hence it follows that if the 
distance from the greatest zero to the point y = 1 is to be a minimum (which must be 
the case for the beam width to be a minimum), then P,_,(y) must be proportional to 
the Chebyshev polynomial T,_,(7). Examples appear in Dolph’s paper. 

Relation of the Chebyshev polynomials to filters. Consider the netavork shown in fig. 6, 


Vo to V e2) In-1 Vn jn 
Xj x; xy Xj 
Xe Xo Zs 








Fic. 6. Finite chain of reactive elements connected in ladder form. 


which is composed of n identical circuits of isolated lossless elements. Let the driving 
emf at the input terminals be a sinusoidal source of amplitude V, and of angular fre- 
quency w. Let x, and x, represent the reactances of the circuit elements at this frequency. 
Then 

%,1, = tif oe rer — 27,) 


|. 9 
T My 


or, by placing x = 2/22 


I, ,— «i, + I, , = 0. (38) 


From equation (38) it can easily be shown that if 


” . iol ] n—2 a 2 n-4 
y,(@) = 2 - ( I ); +( 9 ): — se 


then in terms of J, an J, , J,,; may be expressed as 
Tees = Yel — Yr-rlo - (39) 
In terms of J, and V,: 


V;, = —Yp-1 X11 9 + (Yu-1 = Yr-2) Vo ’ 


I, = (Y, = UP alle —_ Yx-1V 0/22 . 


The y’s are polynomials satisfying the same recurrence relation as (38). By comparing 
the expansion of y,(a) with the expansion (22), the following relation between the y's 
and the Chebyshev polynomials can be derived: 

+ 7,(x/2) if n is odd.] 


y,(x) = ol Tax 2) + T,-2(x/2) + -:> (40) 
+ 7,(x/2) if n is even.] 
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Through the use of equation (40) not only can the successive currents and voltages 
in the network shown in fig. 6 be expressed in terms of Chebyshev polynomials, but also 
the input impedance and transfer impedance for any load z, . 

A pplication to directional couplers. It has been pointed out by F. Bolinder 
one wishes to maximize the ratio of the width of the pass band to the level of the reverse 
wave in the pass band in a directional coupler, then for a given number of elements one 
will choose the coupling factors in such a way that the variation of the reversed current 


(2) that if 


level with distance will be a Chebyshev polynomial. 

Conclusion. This paper has presented the properties of the Chebyshev approxima- 
tion method and the Chebyshev polynomials which may be of interest and use to eng- 
neers and physicists. The application of the method to filter design has not been covered 
here because of the length of the subject and the detail required. Glowatzki (16) gives tables 
for determining the number and sizes of elements needed in a filter which meets certain 
specified requirements. A number of other references on this subject are given in the 
bibliography, where also may be found references to articles of a detailed mathematical 


nature for those interested in this aspect of the problem: 
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THE PLASTIC INSTABILITY OF PLATES* 


BY 
H. G. HOPKINS** 


Manchester University 


1. Introduction. A flow theory of the plastic instability of plates under edge stresses 
which can be used to provide information for structural design is not yet available. In this 
paper, an analysis of this problem, for rectangular plastes under uniform edge direct 
stresses, is given that provides information of a general nature only. In comparison with 
previous work by Handelman and Prager (1), Hopkins (2) and Pearson (3), the present 
paper exhibits both similarities and differences. In the first place, Refs. 1-3 and the present 
paper develop analyses within the framework of the elementary bending theory of plates 
and based upon a plastic flow theory relating to an idealized material whose mechanical 
behaviour is that of a work-hardening material obeying the Mises plasticity condition; 
further, in all these papers, the applied edge direct stress distributions are similar, always 
being uniform along an edge and taken to vanish for one pair of opposite edges in Refs. 
1 and 3. In the second place, there is an implicit neglect of strain compatibility in the 
plate middle surface in Refs. 1-3. In brief, in comparison with the analyses of Refs. 1-3, the 
important features of the present analysis are that strain compatibility of the plate 
middle surface is included and that the development of the stress-rate/strain-rate re- 
lationships is more incisive. The present paper is based upon Ref. 4, and certain parts 
of the analysis developed there are here either curtailed or omitted for conciseness. 

2. Notation. The length, breadth and thickness of the plate are respectively a, 6 
and 2h. The tensor notation employed in Section 3 is described there. In Section 4, the 
origin O is a point common to an edge of length 6 and the middle surface; the z- and 
y-axes are in the middle surface, and parallel to the plate edges; and the z-axis is such 
that the frame of reference O(x, y, z) is right-handed. The displacement of a point 
during instability has the components hu, hv and hw respectively parallel to the z-, 
y-and z-axes. The notation for strains, stresses and plate resultants is standard except 
that reduced stresses, i.e. stresses divided by Young’s modulus E, , are used. The tangent 
modulus of elasticity is FE, and Poisson’s ratio (in the elastic range) is (0 < » < 4). 
The analysis is simplified with the use of the non-dimensional coordinates £, » and ¢ 
defined by 

t= 2/a, n = y/a, ¢ = 2/h, (1) 
and then, with an appropriate choice of the co-ordinate system, the origin now being 
taken at a corner of the plate, the unstrained plate occupies the region defined by 


O<E<1, O<Sn<b/f, -1<S¢<1. (2) 


Immediately before normal displacements occur, the plate is subject to the state of 
uniform normal edge stress defined by 


o, = —o,(¢ = 0,1; 0<7< b/a; —1<¢<)D), 
(3) 
o, = —o,n=0,b/a; OS ES]; -I1s fs 0, 
"Received July 24, 1952. 
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where at least one of co, and a, is necessarily positive if normal displacement of the 
middle surface is to be possible. The following abbreviations are used: 


a=a b. 8B = h*/a 3 


of = 3(1 — v)o,/8 a, = 3(1 — v’)o./B 
r, = (26, —o a, — 2e = (2c, — a,)/(o. — 2a,), | 
p= -v-Di-z alent 4 
\ ? r r; Fi Ye | 
=(1 —» + (A, 1yf{l + 1/7, — 2y/r,} > O, 
Cc, = (A 1)/L(20, — o,) c. = (A. — 1)/L(20. — @,), 
| 
oh F » ) | 
w=1-2/s(1- (A, — 1 1 - AX — I) 
/ ( a +( ~} | 
(4) 
= —] — Al —v)/(A, — 1I)(1 + 1/75 — Wr) < -1 
V = {(1 ri), — 1) + (l —v/r,)(Q2 — 1)}/401 — »’) | 
| 
= (A (1 + 1/7; — 20/7 1/1 —v) > 0 
d* = 1(2 T ( l ? | 
6 ad*(x _ 1)(] = P/F I 6 = —<d*(r ae 1)(1 ~ w/t) / hy | 
Dy, =1+ v/r,)8y | 
D.. = 1+ (1 — v/r.)é 
| 
2),.=2+ 0 — 1/7,)6, + & — 1/r.) 6 . 


In the analysis it is convenient to treat with rates-of-change of certain mathematical 


quantities rather than with their small changes, and primes generally denote time 
differentation. 

3. Relationships Between the Stress- and Strain-Rates. In this Section, following Refs. 
(5) and (6), relationships are determined between stress- and strain-rates for a special type 
of work-hardening material; these relationships are developed within the framework of 
the elementary bending theory of thin plates and are appropriate to a rectangular plate, 
loaded initially only by uniform direct stresses along its edges, and in a configuration 
close to its initially flat configuration. Tensor notation is used in this Section; thus, if 
Ox, denotes a right-handed rectangular cartesian frame of reference, then o,;; denotes 
the stress tensor and e’; denotes the strain-rate tensor corresponding to the stress-rate 
tensor o’,, where Latin suffixes have the values 1, 2 and 3, and the repeated suffix con- 
vention for summation is employed. For simplicity, reduced stresses and stress-rates, 
i.e. stresses and stress-rates divided by Young’s modulus E, , are employed. 
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The plastic strain-rate «/? is defined by the equation 


§= 6-65, (5) 
where the elastic strain-rate ¢/* is related to the stress-rate o/; through Hooke’s law, viz. 
ef = (1 + vot; — vo.d;; « (6) 
Then it is proved in Ref. 6 that the plastic strain-rate is given by the equation 
a GZ Ft (7) 
O6;; OC~%,) 


where G is a non-zero scalar that may depend upon the current stress- and strain-states, 
together with their histories, but does not involve the current stress-rate, and where f 
is the loading function. Thus, from Eqs. (5)-(7), the stress-rate/strain-rate relationships 


are given by the equations, 


F : ; or ee : of 
ef, = (1 + v0); — vod; + CG = = on aw oh > @, 
Og; ; 067%; OC>“K) 
(8) 
; , , of 
e, = (1 + v)0;; — vond;; , =~ @, < ©. 
OC >: 
It is assumed now that the plastic strain-rates are incompressible so that 
e. = 0, (9) 
and from Eq. (7) it follows that 
of of of of 
ae = a es — = () 10 
OO“. Oo; 0022 + 0033 ; 
or, equivalently, 
rf I 
“- =Q@ where C=-on; (11) 
Oa 3 


and hence the loading function is independent of o and therefore depends only upon 


o.* where 


ee ee (12) 


G;; = ¢ 
j ij 3 


If the co-ordinate axes Ox, are suitably oriented then the initial state of stress in the 
plate, i.e. just before normal displacements occur, is represented by the tensor 


* 0 0 


, 0}. (13) 


a o 


0 0 0 


In the present connection, i.e. within the framework of the elementary bending theory 
of thin plates, the stress-rate tensor to which explicit attention is given when normal 


displacements first occur is represented by 


0}. (14) 
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Therefore, from Eqs. (8) and (9), the stress-rate/strain-rate relationships, significant 


in the present connection, are given by the equations, 


OF @f , 


, i | y J 
€ag = (1 + veep — VO,,5ag + Gs a Orxs »| 
O0 ag OC,5 me 
of 
| . 
=r a 0c, 
; af of\ af , a“ 
6 = —VvO G ay TT . 5 
Oo Oo Oo 
€,.2 = T Vid = Vo . OaB . 
Oo} i 
ae <o 
, Oo. 
€ ( 


O43 > 0, | 


(15) 


where Greek suffixes have the values 1 and 2, and the repeated suffix convention for 


summation is employed. It is reasonable to assume that the direct stress-rates produce 


zero or negligible shear strain-rates, and then, from Eqs. (15), 


= tor o = Co =o = @. 


Eqs. (15) now simplify to the equation 


I 


€ 7 Vo 
\ 
Oo} f of P of; > 
a G - 1 " O22], 
Og Oo Oo / 
“74 1 + v)o! 
6. = —ve', +o of of 
—— ¢i, + = es: > 0, 
00711 Od» 
0 oT oO} , | 
a G (- : ee + oO 
Oo OG7*: Oo / 
, , ' , 
€ —V(O11, - Foo) 
Of of Of Oo} 
—~ GA 4 SLY a, + “2 a:,), 
Og O020/ \O0;; O00. 
, , , 
€i) = 011 — vo 
| 
, , , 
eo. =(1+p)¢ | of 
6] Of 
= - + a Soe < 0 
r r , 00}; Oo 
€ = “7 43 oC 
es, = —v(ci, + o,), 


Now define two scalars \, and \, through the equations 





(16) 
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af 
eo er ae ge = 0 and = i. > 0,| 
005; i 
(18) 
af | 
» =e./e, of =0 and ——ei,>0,] 
062 J 
and then, from Iéqs. (17), 
a . 2 a) ° 2 
4, = 1+ ae J ), » = 1+ of 2! y, (19) 
00), 0022 
and hence 
oy, ; 9 
Gg? 2 eaQg,-2”, G®~ = +0, - 1)”, (20) 
0011 0022 


where, in the first equation, the positive and negative signs correspond respectively to 
loading and unloading in the special case o{, > 0 with o3, = 0, and similarly for the 
second of these equations. Therefore, from Eqs. (17) and (20), it is possible to formulate 
stress-rate /strain-rate relationships that involve the two scalars \, and \, whose numerical 
values could be determined from experimental data. Note that 


af ar 2 
(Ay — 1)/Q2 — 1) = (2f / aL. - 


The simplest possible assumption for f is that the Mises plasticity condition applies, 


V1Z 


(22) 


where the strain-history enters only through the single parameter k. Loading occurs 


only if 
f’ = ot,o;* > 0, (23) 
and in the present connection the criterion is therefore that 


f \f 
\O1u1 Ta an (o1, -+- 0G 2) f' oF) ‘oa 
\ o ‘ e 


(oi, + zy 


( ) 
} l \ 1 , | 
+ )o22 — 3 (on + oud) et — 3 (on + xia} + 9 (on + o22)(oi1 + 022) > 0, (24) 


\ 
i.e. 
; I ; ; 
id = 3 {(20,, ad Ceel@is + (2622 —_ 011) 042} > 0. (25) 
Thus 
Of _ 1, of _ 1 
de. = 3 (26, O22), OGes = 3 (2¢22 = 1), (26) 
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and accordingly, from Eq. (21), 
(A, — D/As — I) = {(201 — o22)/(2022 — o1,)}’. (27) 
From Eqs. (19) and (26), 


G = 9A, — 1)/(2011 — O22)? = WA2 — 1)/(2o22 — O11) 5 (23) 


and G is calculable if either \, or \. is known from experimental data. It is possible 
therefore to formulate the stress-rate/strain-rate relationships in terms of a single 
scalar to be determined necessarily from experimental data; these relationships are 





; a ty 
é€f, = Aon — jy + p> (Ai — 1) fo22 
\ R, } 
éio = (1 + vj)oi2 , 
io = —lv + > (Ae — I) fons + Aroine , ‘ or | 
: 0) y) (261; — O22)o011 + (2022 — 711)022 > O, 
, - is, 
j= —jy tA. —1 —Z (Ce — Yee 
R, 
r (29) 
i 4 2, 
—\v+ (A. — 1) — BF A: — 1) /o22 
R, 
— , , | 
oy C1, ~~ Pen 
} 
eo = (1+ »)012, | 
f 2611 — F22)01; + (2022 — a1) 022 < 9, 
62 = —voi1 + O22, | 
| 
| 
63 = —v(oi, + o22),J 
where 
R, = (261: — 622)/(611 — 2622), Rz = (2022 — o11:)/(¢22 — 201:), (30) 


and , and 2, are related through Eq. (27). These stress-rate/strain-rate relationships 
agree with those obtained in Refs. 1 and 2. However, in Ref. 2, the result Eq. (27) could 
be deduced only for the particular case discussed in Ref. 1, viz. o22 equal to zero, and 


the result is then 


r 


> 


] ‘ ‘ 
= 4 (hi + 3), (31) 


where X, is interpreted simply as E,/E. 
In the next Section, it is more convenient to use an engineering notation for stresses 
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and strains, and accordingly Eqs. (29) are rewritten now in the form, 


l 
ef = \,0. — jv + — a = I) fou, 








vt, = U1 +»)z 
if ex aa ee — 1) el + dso? | 
; Ms f . saci (20, 7 o,)0% + (20, = a,)o, < 0, 
f L 
6 = -\Vv +r (Xr, = 1) == (Ay — 1) ?o; 
Ns J 
f (32) 
1 ! 
—<y t+ (rh, — 1) -—-- — I) fos 
ry 
é’ a’ — vo 
vi, = Al + »)71,, 
((20, — o,)o; + (20, — a,)o, > 0. 
e = —volL +a’, | 
e = —r(ei + o),) J 


Let the rate at which work is done by existing stresses on the change of shape of an 
element be |W’ per unit volume so that 


W’ = E,0;;e'3% . (33) 
It may be proved from Eqs. (8) and (22) that 
((1 +») + 2Gk)Bf’, f'> 0,| 


W’ = 
| (1 + Ef’, f’ <0.) 


(34) 


Thus the loading condition that /’ is positive is interpreted that W’ is positive; the 
neutral condition that /’ is zero is interpreted that W’ is zero; and the unloading condition 
that /’ is negative is interpreted that W’ is negative. From Eqs. (13) and (32) the load- 
ing and unloading criteria are now written in the form 


{(2o, — o,)v + (20; — o2)}e2 + {(20, — o.)v + (20, — a) }eh $ 0, } 
(35) 


Ww’ >0, 


respectively. 

4. The Fundamental Equations. In this Section, the fundamental equations are de- 
veloped on the basis of the stress-rate/strain-rate relationships in Eqs. (32). It is assumed, 
exactly as in elastic thin plate theory, that plate elements originally normal to the middle 
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surface remain normal to the displaced middle surface. Then 


eo = ef — Brew’; , | 
| 
| 

6, = « — Biv’,,, (36) 
; ‘ ™ | 
V1 = 7’ — 2Bew’,, ,) 


where ¢/ , «s and y’ are strain-rates in the middle surface. Here, and also elsewhere in 
this paper, suffixes ¢ and 7 (and also n and s) preceded by a comma denote partial differ- 
ential coefficients. From condition (35) and Eqs. (36), the loading criterion is 


ft 2} fo if K’ 20 (37) 
respectively, where 

Bek’ = (26, — o2){(1 — v/ref + (v — 1/r,)e}, ) 
‘ (38) 

K’ = (20, — o2){(1 — v/r,)w’es + (Vv — 1 r;)w’,»}.) 
It is necessary always that — 1 < { < 1;if — 1 < ¢ < 1, then the surface with equa- 
tion ¢ = ¢, separates loading and unloading regions, and on this surface W’ = 0; f) = 
+ 1 corresponds to either loading or unloading over the entire plate thickness, and on 
the (plate) surface ¢ = + 1, in general W’ + 0, and therefore this surface is not in this 


sense a boundary to a loading or an unloading region. It is seen by setting either w’,, 
or w’,, equal to zero that the sign of K’ is an indication of the sense of bending. 
The resultant force-rates N/ , Ni and N12, are defined by the equations 
SRT r ry , , , - : 
(N2,N2,N1)/h= | (02,02, 71, dt, (39) 


and therefore 


al 


(Nt —vNZ,Ni—wND/h= | (0! — vol, of — vo!) dt. (40) 


/-1 


From Eqs. (32) and (38) and conditions (37), 
; ions * 
of — vo, = € + Beil(E — fo) K’, 





of — vor = €, + Bef — fo) K", (W’ > 0), (41) 
ri, = ¥!,/2(1 + »), } 
and 
oc —vl=eé 
| 
of —vwl=e,, ( (W’ <0). (42) 


rl, = v!,/21 + »),) 
Conditions (37) show that there are two cases to be considered and the analysis now 
proceeds differently according as K’ > 0 (Case A) or K’ < 0 (Case B). 
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Case A. Here K’ > 0, and then (1 > ¢ > &) and (¢ > ¢ > — 1) are respectively 
the loading and unloading regions. 
From Eqs. (36) and (39)-(42), 


(NL — wN1)/2h = & + 46c,(1 — §)’K’, 


(Ni — vN!)/2h = & + 38c.(1 — &)’K’, (43) 


Ni,/2h = y'/21 +», 
and therefore 
—38(1 — f0)*K’ = fet — (NI — wN{)/2h}/e, = [e — (Ni — wN2)/2h}/e.. (44) 
From previous remarks, and from Eqs. (38.1) and (44), f% = ¢ if — 1< ¢ < 1 and 


¢, = + 1 otherwise where ¢; = ¢*(P) is defined by the equation 

c* — 2Me* +P =0, (45) 
i. 

{* = M2 (M’* — P)™, (46) 


and 


P = 1 — 4L{(20, — o.)(1 — v/r,)pi + (202 — o,)(1 — v/r2)po} 


+ {( — 1)(1 — »/r,) + A, — DL — v/r)}, 
r (47) 
p, = (N! — »N‘)/2hpK’, 





D2 = (Ni — wN1)/2hBK’. 
If — 1 < ¢*(P) < 1, where ¢*(P) is defined by Eq. (46), then ¢, = ¢*(P) and the surface 
with equation ¢ = {) = ¢*(P) separates loading and unloading regions of the plate. If 
such a surface of separation exists in the immediate neighbourhood of some particular 
point of the middle surface, then, at this point, it is necessary first that 


M’?>P (48) 
and second that it is then possible to choose the sign in Eq. (46) so that 
~i Ss fh. Ss 4; (49) 


if either of the conditions (48) and (49) is violated, then loading or unloading occurs 
over the whole of the plate thickness in the immediate neighbourhood of this point, and 
then ¢) = + 1. It is therefore not impossible for loading (or unloading) to occur over 
the entire thickness of part of the plate and a varying degree of loading (or unloading) 
to occur over the remaining part. Note that both roots of Eq. (45) may satisfy — 1 < 
¢*(P) < 1 in which case there are two possible solutions. 

The resultant bending and twisting moment-rates M! , Mj and M!,(= — M},) are 
defined by the equations 


1 
(Mi, Mi, Mt)/h? = [ (ot, of, — rhe as, (50) 
1 
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and therefore 


(Mi — »M! , Mi — »M!)/h = | (of — vol, of — vo'l)e df. (51) 


From Eqs. (36), (41), (42), (50) and (51), 


Mi — »vMj = —§h'B(w’s: — de, K’), 
Mi — pM! = —2h'B(w’,, — de.K’), (52) 
M!, = —Mi, = {3h’B/(1 — v)}(1 — vw’, , 


where d = d*(f,), and therefore 


—M2/{Rh'B/(1 — v*)} = {1 + 61 — v/r) jules + tv + 600 — 1/12) } Way || 
| 


—M!/{2h*B/(1 — v)} = fv + 6 — 1/7) u's + {1 + 6.01 — v/rr) fw, 5 (53) 
M1,/{3h?B/(1 — »*)} = (1 — vw" - 


The resultant total vertical shear force-rates V! and V/ are defined by the equations 


aV! = M!, — 2M!,,, | 


(54) 
aV; = M,, — 2M,..,) 
and therefore from Eqs. (53) 
— V!/{2ap’/(1 — v’)} =5(0 + 6,(1 — v/r,)}w’:-+{(2— v) + 6,(v — 1/r.)}u’,,], 
(55) 
— Vi/{2ap°/(1 — v)} =5et1 + 6(1 — v/r.)}w’,,+{(2 — vy) + 6, — 1/7) } uw’). 


Now consider the equilibrium of the plate when small normal displacements have 


§ 57) of the plate, for small normal displace- 


occurred. The equations of equilibrium (7, § 
ments of its middle surface, are 


Mrs — 2Miyeg + My igy = —h(NW,2¢ + 2N2W.t_ + Nias); 
Noe + Nz... = 9, 4 (56) 


Nye + Ny, = 0; 
partial differentiation of Eqs. (56) with respect to time shows that 


M! we — 2M!) ey +H Mi gn = 2h%(o,w! ge + o2W' a), 
Nie + Niu. = 9, (57) 


Nie +N’, =0, 


ry,t 


immediately following the incidence of normal displacements because immediately 
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before these occur w and N,, are identically zero. From Eqs. (53) and (57.1), 


oe 9° a” 
V iw + (1 — »/7) = (6,w’ we) + — T/rs) 53 (620% 


oF dg 
‘ee a a” nee . o- vr _ Ks ron? rn 
+ (v 1/r:) —s (6,w’es) + (1 — v/re) <3 (6,0) = —(otwlee + o30',,), (58) 
On - On ‘ 
where 
0” AY 
Viz ( ; :) ; 5f 
: ry + op (59) 


Equations (57.2) and (57.3) are satisfied identically through the introduction of a 
stress-rate function ¢g’ that is related to the resultant force-rates by the equations, 


N2/2h = o' mn 5 Ni/2h = ¢' et ; N!,/2h = —o':, - (60) 
The strain-rates ¢! , ef and y!, are given by the equations 
= B7w,, &§ =B',, w= B°w, +H), (61) 
and are compatible only if 
, > 
Exon + €) te i Vouste = 0, (62) 
i.e. from Eqs. (36), 
Ge + €2,8¢ 7 V.20 oe 0, (63) 


which is an equation that involves the middle surface strain-rates only. From Eqs. (44) 
and (60), 


€{ a Pian - ve" ee sie 7c, (1 = fo) K’, 
2 = lee — Way — F8c(1 — f0)K’, (64) 
. y= -Al + ele, 
and therefore Eq. (63) 1s 
4 Ff l 3” 2u-7 
Vie’ = 78 ao ta —3}{(1 — f)°K’}. (65) 


Case B. Here K’ < 0, and then (1 > ¢ > g) and ({& > ¢ > — 1) are respectively 
the loading and unloading regions. It is straightforward to show that the analysis of 
Case B is deducible from the analysis of Case A if all auantities are now referred to the 
coordinate system O’(t’, n’, ¢’) where 


£’ = const. — &, n = const. — 9, i= -f, (66) 


the values of the constants being irrelevant. Full details are given in Ref. 4, and here 
only the results are stated. The equations that correspond to Eqs. (53), (55) and (58) 
are identical in form but now d = d*(— ¢,) and ¢, = — ¢*(— p, , — Po); the equations 
that correspond to Eqs. (60) are identical in form; and the equation that corresponds to 
Eq. (65) is 


a 


hse l 
Ve =- ac. & eto 5e 


; JA + °K") (67) 
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where again ¢, = — ¢*(—p, , —p.). The convention is now adopted that the above 
changes are implied if A’ < 0, and then Eqs. (53), (55,) (58), (60) and (65) apply un- 
restrictedly. 

The fundamental mathematical problem is therefore to solve the non-linear, de- 
pendent, fourth order, homogeneous, partial differential equations (58) and (65) subject 
to boundary conditions on w’ and g’. These boundary conditions depend upon respec- 
tively the plate edge support and the applied resultant force-rates, and, in general, are 
not arbitrarily prescribed. If the middle surface is sub-divided into two or more regions 
by a curve c, say with normal direction n, on which K’ = 0, then it is necessary on 
physical grounds that w’ and g’, together with their first three partial differential co- 
efficients with respect to n, are continuous on c; but, from Eqs. (58) and (60), Ww’ nnn 
and @’nnnn are not, in general, continuous on c and therefore in general the analytical 
form of the solution is not the same on either side of c. 

Now consider the application of the analysis to two particular classes of problem. 
These are often designated in the published literature as the von Karman and Shanley 
problems; in this paper they are designated briefly as problems A and B, respectively. 

Let n and s denote directions respectively normal to, and along, a plate edge. 

Problem A. Here it is postulated that the changes of the applied edge stresses at 
instability are zero. Then, by definition of the problem, 

Gs = Ons = 0 (68) 
along an edge, but note that, because Eq. (65) is not homogeneous in gy’, Eqs. (68) do 
not imply that ¢’ can be taken as identically zero. The boundary conditions on w’ depend 
upon the plate edge support; for example, if the edges are either clamped or simply- 
supported (which are important practical cases) then the boundary conditions on w’ 
are respectively either 


w= Ww =: () (69) 
or, from Eqs. (53), 
w’ = w’,, = 0 


along an edge. The solution of Eqs. (58) and (65) subject to the boundary conditions 
(68) and, for example, (69) is clearly, in general, an extremely difficult problem. If the 
instability stresses do not greatly exceed the elastic limit, then it is not unreasonable 
to approximate Eq. (65) by the equation 

Vie’ a= (). (70) 
This approximation is made implicitly in Refs. 1 and 2, and results in a considerable 
simplification of the problem. The appropriate solution of Eq. (70) subject to the con- 
ditions (68) is 

gy’ = 0. (71) 
Then, from Eqs. (46) and (47), remembering that JM is not greater than — 1, 

f= M+ (M — 1)” (72) 


and this equation necessarily gives a value for ¢* such that — 1 < ¢* < 1 and therefore 
tf. = + ¢* according as K’ 2 0; there is now a surface separating loading and unloading 
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regions and it has the equation ¢ = + ¢* according as K’ 2 0, i.e. a surface composed 
of parts of the planes ¢ = + ¢* with discontinuities at points such that K’ = 0. Such 
discontinuities may be due to the approximate nature of the solution. Further d = 
d*(+¢,) according as K’ 2 0 and therefore d = d*(¢*) in both cases. Equation (58) is 


how 


Dw eeee + 2DpqW' tay + DoW ayn. = —(oiwiee + o3W' yy) (73) 


where D,, , Doo and D,, (see Eqs. (4)) are constants dependent upon the stress-state. 
The problem therefore depends upon the solution of the linear, fourth order, homo- 
geneous, partial differential equation (73) with constant coefficients dependent upon 
the stress-state, subject to boundary conditions on w’, e.g. Eqs. (69); the characteristic- 
values and characteristic-functions of this equation correspond respectively to the 
instability stresses and normal displacement-rates. Note that, although Eq. (73) has a 
doubly-infinite set of characteristic-values (and corresponding characteristic-functions), 
only the first is of importance from a practical standpoint, and further only in this case 
may Eq. (70) be a useful approximation to Eq. (65). Equation (73) is given in Ref. 2 
and also, for the particular case o, equal to zero, in Ref. 2. As remarked in Ref. 2, there 
is no serious mathematical difficulty involved in the determination of the (approximate) 
instability stresses for the edge conditions that are of practical importance; indeed the 
only difficulties are in the computation, and it is useful to note that it is simplest to 
regard o, , o, and @ (or 8) as known and to solve for 8 (or a). 

Problem B. Here it is postulated that the changes of the applied edge stresses at 
instability are such that plastic deformation only occurs. 

It is evident on physical grounds that, if it is possible for loading, simulateously 
with bending, to occur over the whole of the plate, then the plate bending stiffness is 
least and therefore so also are the instability stresses. An analysis based on this assump- 
tion, in any case, must lead necessarily to a lower bound for instability stresses; the 
assumption is discussed later. 

Then, by definition of the problem, {, = * 1 according as K’ 2 0; notice, however, 
that the (plate) surface £ = + 1, in general, is not a boundary to a loading or an un- 
loading region. Further d = d*(— 1) = 1 in both cases. Equation (58) is now 


Dy, eeee + 2Di2W' geny + DoW’ anny = —(oiw’ ee + o2W' yy) (74) 


where D,, , D.. and D,, (given from Eqs. (4) after setting d = 1) are constants dependent 
upon the stress state. Equations (73) and (74) are identical in form, and therefore the 
determination of the instability stresses in this problem is a task exactly comparable 
with that of the determination of the instability stresses in the previous problem. Equa- 
tion (74) is given in Ref. 2 and also, for the particular case o, equal to zero, in Ref. 3. 

In the discussion of the assumption that plastic deformation only occurs, the analysis 
proceeds differently according as K’ > 0 (Case A) or K’ < 0 (Case B). 

Case A. From Eq. (65) it follows that ¢’ satisfies the equation 


Vie’ = Bi Sry w' eee rs (6, + 52) W' tenn + Sor2 'W' sane} (75) 


and satisfies boundary conditions necessarily consistent with the assumption that load- 
ing occurs everywhere. If this assumption is true, then at every point of the middle 
surface, it is untrue that — 1 < ¢*(P) < 1. The critical case occurs when {*(P) = — 1. 
Now suppose that, in general, at some particular point of the middle surface ¢*(P) = 


























198 H. G. HOPKINS {Vol. XI, No. 2 


—(1 + 5); for simplicity, it is assumed that 6 is real. Then, from Eq. (45), it is straight- 


forward to prove that 


N& —6+6,=0 (76) 
where 
: j : . iy K’ — 
oy = The: — 20,)¢ te + (G2 — 2a:)¢ nn) re} a — ] (44) 
and hence 6 > 6, as N > O. The critical case will occur if 6 = 0 and then 6, = 0. A 


sufficient condition for loading to occur over the entire thickness of the plate is that 


5, > 0, i.e. from Eq. (77) that 


p ty e 
H(o, — 2or)giez: + (o2 — 2a en} / (8 i ) oe & (73) 


Case B. It is straightforward to show that the analysis of Case B is deducible from 
the analysis of Case A. Full details are given in Ref. 4, and here only the results are 


stated. The relation corresponding to relation (78) is 


F K’ a 
{(o, — 2o2)¢’r: + (a2 — Zoey wt /(-6 A) > (79) 


Therefore a sufficient condition, applicable in either case, for loading to occur over 


the entire thickness of the plate is that 
oe , ee es 
{(o, — 2o.)e’es + (o2 — 201) Gan} L K ae OE (80) 


Here it is required that loading occurs over the whole of the plate. Therefore it sufficient 
that ¢’ is determined from Eq. (75) subject to boundary conditions consistent with the 
condition (80) for loading everywhere. The accurate determination of yg’ may not always 
be very simple. However, the following approximate analysis indicates that it should 
always be possible to choose the rate-of-change of applied edge stress at instability so 
that loading occurs over the entire plate, although the determination of the least possible 
rates of change is, in general, a difficult problem. Now if the instability stresses do not 
greatly exceed the elastic limit, then it is not unreasonable to approximate Eq. (75) by 
the equation 


Vie’ = 0. (81) 


Next take the boundary conditions on ¢’ to be that (i) ¢’,, is constant along any edge 
and has the same value for opposite edges and (ii) ¢’,, is zero along all edges. The solu- 
tion of Eq. (81) subject to these boundary conditions corresponds to a uniform stress- 
rate distribution. The condition (80) is now simply 


const. /8 oS 2 (82) 


and loading will occur over the entire plate thickness if of = — g’,, and oj = — g’¢; 


are such that 


ie eet +, obo F Max. | K’ |. (83) 

















1953 THI PLASTIC INSTABILITY OF PLATES 199 


The above analysis of Problems A and B refers to two salient problems of the general 
analysis; the first corresponds to greatest possible instability stresses and the second 
corresponds to least possible instability stresses; and, between these extremes, it appears 
probable that there exists a continuous set of instability stresses, this corresponding to a 
continuous set of changes of the applied edge stresses. 

In Refs. 1, 2 and 3 some examples of Problems A, A and B, and B, respectively are 
considered, but in all cases there is implicit neglect of strain compatibility. In Ref. 4, 
a detailed analysis is given of two cases, viz. (i) an infinitely-long simply-supported 
plate under uniform compression stress along one pair of edges and (ii) a finite simply- 
supported plate under uniform compression stress along one pair of opposite edges. 
Problems A and B are solved first accurately and then approximately (through neglect 
of strain compatibility) in case (i) and are solved approximately in case (ii). The approxi- 
mate solutions may be expected to approximate closely the accurate solutions when the 
instability stresses are not too greatly in excess of the elastic limit. 

5. Conclusions. There are two salient problems, often designated in the published 
literature as the von Karman and Shanley problems; in the first, which corresponds to 
greatest possible instability stresses, it is postulated that the changes in the applied 
edge stresses at instability are zero; and in the second, which corresponds to least possible 
instability stresses, it is postulated that the changes in the applied edge stresses at 
instability are such that plastic deformation only occurs. It appears probable that, 
between these extremes, there exists a continuous set of instability stresses, this corres- 
ponding to a continuous set of changes in the applied edge stresses. 

In the first problem, if strain compatibility is neglected, then the determination of 
the instability stresses is straightforward; it is probable that these approximate values 
underestimate the accurate values on the basis of the present theory by an amount that 
is small only if the instability stresses do not greatly exceed the elastic limit. In the 
second problem, the determination of the instability stresses is straightforward; however, 
although the assumption that only plastic deformation occurs at instability is probably 
correct, the determination of the corresponding least possible changes in the applied 
edge stresses is in general difficult. 

It is known that theorectical values, based on the present flow theory, considerably 
overestimate experimental values of plastic buckling loads for long simply-supported 
plates under end compression. Further, deformation theories are available that provide 
results in good agreement with the experimental results; however, this fact is sometimes 
taken as a justification for a deformation, as opposed to a flow, theory of plasticity, an 
such a justification is untenable. If the experimental results relate to perfect test speci- 
mens perfectly loaded or, if these results are not particularly sensitive to imperfections 
of geometry, loading or mechanical behaviour, then the disagreement should be trace- 
able to the assumed mechanical behaviour; otherwise, the disagreement should be traced 
to either imperfections or assumed mechanical behaviour or both of these matters. 
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WATER WAVES OVER A CHANNEL OF INFINITE DEPTH? 


By 
THOM R. GREENE (The General Electric Co., Schenectady, N. Y.) 
AnD ALBERT E. HEINS (Carnegie Institute of Technology) 


I. Introduction. This is a continuation of some work on gravity waves in channels 
of finite depth with special types of obstacles (1), (2), (3). Thus far three cases have 
been considered: the channel of finite depth with (i) a dock; (ii) a submerged plane 
barrier and (iii) the joining of two free surfaces. In this paper and one to follow we propose 
to discuss the corresponding problems for channels of infinite depth. The problem which 
concerns us here is case (ii) and its limiting case (i) for such channels. 

The physical background for the problem has been discussed in (1) and (4). Here we 
shall reformulate the problem so that the notion of infinite depth is admitted. As we 
have done in the past, only propagating regions will be considered here, the possibility 
of non-propagating regions may be treated in a fashion similar to that of the propagating 
region. It will be assumed that the waves have propagation normals which are not 
perpendicular to the edge of the barrier and in the beginning we shall assume that the 
barrier is submerged (2). As limiting cases, we will get the case of normal incidence and 
the ‘dock problem.”’ 

While there is a certain parallelism between the formulation of the problem in channels 
of infinite and finite depth, we observe that the mathematical methods of solution become 
more subtle in the channels of infinite depth. Still, we provide here methods which will 
solve the problem we have described. With minor modifications in method but definite 
changes in formulation, we can solve (iii) for infinite channels. 

II. Fundamental Equations of Motion and Statement of Problem. If we use the linear- 
ized equations of motion as we have done in the past, we have for an incompressible, 
non-viscous fluid, in irrotational motion that 


o:2 + dy +o. = 0 
where ¢’ is the velocity potential. On a rigid surface ¢, = 0, while on a “free surface” 
gy, + dit = 0. 
Assuming monchromatic time dependence 
¢’ = $(x, y, z) exp (—iwl) 


and these conditions become ¢, = 0 and g¢, = w'¢ respectively. 

Here we are concerned with a channel of infinite depth with a semi-infinite plane 
barrier submerged a units below the undisturbed free surface (See Figure 1). We shall 
assume that wave motion exists for | x | —~«, y = 0 where the propagation normal is 
parallel to the free surface but not at right angles to the edge of the barrier, that is the 
2 axis. As such there will be a z variation exp(ikz) (k > 0) in @ which we may suppress 
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from the beginning. This gives us as fundamental equations of motion 


¢,,+ ¢, —ko=0 in the fluid medium, (2.1) 
, = 0 on a rigid surface, and (3.1") 
o, = Bo B =w'/qg on a free surface. (23/5 


Equation (2.1) possesses the following types of bounded solutions. For z >, y < a 
@ = exp [tkz + By + 7{8 — hk} a] 
while forz ~—-0~,0O<y<a 
exp (tkz + ixx) cosh ayy/a 

where 

x = {(a,/a) — hk? }'” 
and 

a sinh a = Ba cosh a (2.2) 


has the smallest real roots + a . Eq. (2.2) also has an infinite sequence of imaginary 
roots 7a, ,@-, = —a,,n = 1,2, --- . Wave motion exists if 8 > k and | ap/a| > k. 

III. Derivation of the Integral Equation. Following methods which have been dis- 
cussed elsewhere (1), (2), (3), it is a fairly straightforward task to formulate the integral 
equation which describes the boundary value problem considered in Sec. II. In view of 
the fact that we have suppressed the z dependence in ¢, we have essentially a two dimen- 
sional problem. We shall divide the half space y < a into two regions: (a) the strip 0 < 
y <a, —2© < x <©@ which we shall henceforth refer to as R, and (6) the half plane 
y <0, —2 < x2 <@ which we shall refer to as #, . Since we require the use of the 
Helmholtz theorem in its improper form to provide representations of ¢(z, y) in R, and 
R, in terms of ¢,(z, 0), some statements regarding the growth of ¢ as r © will be in 
order (r = [x” + y’]’’’). To begin with, in R, for x >, (2, y) is asymptotic to 


[y, exp (ixz) + 2 exp (—ixx)] cosh aoy/a. (3.1) 


This is tantamount to the statement that the asymptotic form of ¢ for z —— = is the 
bounded solution which would have been obtained had the barrier extended to infinity 
on the right. For > and y < a, we specify the asymptotic form of ¢ as 


lys exp (tox) + ys exp (—iox)] exp (6y) (3.2) 
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where « = {6° — k’}'” and the positive root is to be understood. This again is a state- 
ment that for a channel of infinite depth without a barrier (3.2) is the bounded solution. 
One final task remains and that is to specify the nature of ¢ asr ~@, 2 < Oandy < 0. 
In this we are not as fortunate as we were in the other two cases. We shall assume (and 
verify with the solution of the problem) that ¢ = 0 [exp (k, r)] for y < 0,7 < 0, r—@ 
where k, < k. This condition will enable us to write a representation for ¢ in R, , that is 
for y < 0. One final assumption is required and that is the statement that ¢ and ¢, are 
integrable for all finite x — again subject to final verification. 

If we introduce a Green’s function G“” (2, y, 2’, y’) associated with Eq. (2.1) and 
which satisfies the boundary conditions 


Gy” = pG"” y = 4, —-oe <xr< o@ (3.3) 


and 
G,’ =0 y = 0, —-e <4< @ (3.3’) 


and the further conditions at infinity that G"’ — 0 for 2’ — x © while for z — 2’ © 


ax( Qa 4 ‘sinh 2a) 


Qn day cosh aoy/a cosh acy’ /a sin K(x — 2x’) 
7 = a — ee ey “¢ 7 


we have as a representation for ¢ in R, 
d(x, y) = [cosh apy/a][y: exp (xx) + y2 exp (—ixz)] 


- | G(x, y, x’, 0),-(2’, 0) dx’. (3.4) 











Here (3) 
- — =. 2a, (cos a,y/a)(cos a,y’/a) exp [—| x — 2’ | {k° + a&/a?}"”7] 
G(2, y, 2’,y) = = a ee es ere 
me a(2a, + sin 2a,){a,/a +k} 
+0 if z’>s or 
4ao(cosh a,y/a)(cosh acy’ /a) sin K(x’ — 2) 7 
—- = : ; f ' , 
+ ax(2ao + sinh 2a) ; ilies 
For the region R, we have in a similar fashion 
o(z, y) = | G(x, y, x’, O)b,-(x’, 0) dx’ (3.5) 


“0 


where 


G""(z, y, 2’, y') = - [Ko(kf(a — 2’ + (y — y’)?}'”) 


+ K(k { (x a x’)? + (y + y’)}'”)]. 
The function K,(4) is the solution of the differential equation 
Ys tys/i-—y =0 


which is O(log 5), 6 > 0°. It is a Bessel function of the second kind which possesses the 
property that it is of the order exp(— 6)/é'” for 6 +. This last order condition ac- 
counts for the choice of the constant k, . On the basis of these two representations of 
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¢(x, y) in R, and R, we can form an integral equation of the Wiener-Hopf type. We note 
that @(z, y) is continuous at y = 0 for all z > 0. We have then for z > 0 


¥1 exp (txx) + y2 exp (—72x2) 
— | [G(a, 0, x’, 0) + G'(2, 0, 2’, 0)]6,(2’, O) dx’ = 0. = (3.6) 


IV. The Fourier Transform of the Integral Equation. We now define Eq. (3.6) to read 
for all x 


I(x) = J(x) - | IG‘? (x, 0, x’, 0) + G(x, 0, x’, O))L(2x’) dz’. (4.1) 
Here 
I(x) = 0 for Se 
J(x) = y, exp (tx) + y2 exp (—7k2) z>v0 
= 0 r <0 
L(x) = 0 z<9@ 
= ¢,-(z) > 9 
Clearly then J(x) = 0 [exp(kz)] for x —~— « and L(x) = O[expk,z]forz—@. Actually 


L(x) is bounded for z —©, being of the form exp(+ icx) but this last description of 
¢(x, 0) for x © need not be stated explicitly. 

Let us now introduce the Fourier transforms of (x), J(x) and L(x) as well as those 
of G (zx, 0, 0, 0) and G™ (z, 0, 0, 0). Henceforth we shall use the notation J(w) to 
denote the Fourier transform of I(x), J(w) for J(x), etc. We have from Eq. (4.1) 


; ¥1 Y2 v(1 (2) 
Iw) = — 4% _ yw) + 6 (w (4.2 
” uw — k) uw + xk) + (w)] ) 
where 
G'"’(w) = | exp (—iwa)G''’(zx, y, 0, y’) dx 
cosh ps[8 sinh p(t — a) + p cosh p(t — a)] 
- pip sinh pa — 6 cosh pa} 
and (s, t) = (y, y’) ify < y’ or (y’, y) ify > y’. Furthermore 
GCG’ (w) = | exp (—iwz)G'”’(z, y, 0, y’) dx 
= fexp[-p|y—y’ |] + exp[—p |y + 9’ |]}/2p. 
where p = {k’ + w’}’”’. The branch of p has been chosen which reduces to k when w 
= 0. 


We now examine the regions of regularity of the various transforms. First /(w) is 
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regular in the upper half plane v > — k (w = u + iv), J(w) is regular in the lower half 
plane v < 0, and L(w) in the lower half plane v < —k, . This follows from the assumed 
integrability for finite x of J(x), L(x) and J(z) as well as their growths for | x | >. 
Similarly G'” (w) is regular in the strip — {k* + aj/a’}'? < v < 0 while G® (w) is 
regular in the strip — k < v < k. There is then a common strip of analyticity for all of 
these transforms, the application of the Fourier transform to Eq. (3.4) is permissible 
and we obtain Eq. (4.2). Upon simplifying we get 


I(w) = y 12 — K(w)L(w) (4.3) 


aw — x) ¥ (w+ x) 
where 
K(w) = |1 — B/p|lexp ap]/[p sinh pa — 8 cosh pa}. 
Following Wiener and Hopf (5), Eq. (4.3) may be rewritten 


K.(w) — K.(«) K.(w) — K(—x) 
K .(w)I(w) - ¥ — = — 
uw + kK) 


zw — K 


= vik (a) + 72K (—«) — L(w)K_(w) (4.4) 
(wo — «) uw + x) ' 
where A ,(w) is regular in the upper half plane v > — k while K_(w) is regular in the lower 
half plane v < 0. Eq. (4.4), asa result of an argument of analytic continuation described 
by Wiener and Hopf, implies that each side if individually an entire function. Our next 
goal is to determine the explicit forms of A_(w) and K,(w) as well as the entire function 
of separation. 

V. The Determination of K_(w) and K,(w). Thus far we have pursued a familiar 
course to arrive at Eq. (4.4). From this point on we encounter a number of unfamiliar 
situations. In the decomposition of A(w) we have three individual components to con- 
sider: 


(a) psinh pa — 8 cosh pa 

(b) exp (pa) 

(c) 1 — B/p 
Component (a) is familiar and has been discussed in (2). The type of factoring we employ 
implies that the portion of (a) which is regular in some lower half plane also has a regular 
reciprocal in the same half plane. The same remark may be addressed to the upper half 
plane factor. We have then 


psinh pa — 6 cosh pa = M,(w)/M_(w) 


where M ,(w) is regular in the upper half plane v > — a,{1 + a’k’/ai}'” and M_(w) 
is regular in the lower half plane v < 0. Here 


M.(w) = [] [{1 + a?k?/az}'” — iaw/a,,] exp (iaw/mm) 
m=1 


and 
= i i cia rae pee oe « > 
1/M_(w) = [] [{1 + a?k?/az}'” + iaw/a,] exp (—iaw/mn) —- (w* — x’). 


m=1 ao 
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For 
w|—> © Imw > —{ai + ak’ }'” 


M .(w) = O[wI'(—iaw/x) exp (—yiaw/m) | 
where y is the Euler-Mascheroni constant, while for | w|—“, Imw < 0 
M_(w) = Ofexp (+7tawy/x)T(iaw/r)w *). 


Now we turn to (b). The decomposition of exp(pa) multiplicatively is equivalent to 
the decomposition of additively. p = {w* + k°}'” is regular in the entire w plane pro- 
vided cuts are introduced along the imaginary axis from k to infinity and — k to negative 
infinity. We define 


1.21 1/2 2, 132 11/2 “16 . 16 
[w+ k°} =|w+k [exp {7/2 arg (w — tk) + 72/2 arg (w + i 
where 


—3nr/2 < arg (w — ik) < 7/2 


—7/2 < arg (w + tk) < 37/2. 


With these definitions at hand we turn to the decomposition of 1/p from which we can 
readily obtain that of p. From Cauchy’s theorem we have 


. 


2ri/(w’ + k’)'? = | dz/(z — w)(z’ + k’)'”” (5.1) 


7C 


where C is a closed rectangular path in the strip — k < Imz < k whose vertical con- 
tributions disappear as they recede to infinity. The upper side of the rectangle may be 
deformed into a path along the branch cut from ¢ to infinity, while the lower one goes 
into a path over the lower cut and we get 


»@ 


wi/(w’ + k’)'? = | dy/(w — iy)(y’ — k*)'? - | dy/(w + iy)(y’® — k’) 

“k “« 
This provides a separation of the required type—the first term being regular in the lower 
half plane v < k, the second in the upper half plane v > — k. For our purposes, however, 
we require more explicit information, particularly as to the nature of these integrals as 
w|—o and we therefore proceed with the evaluation of these integrals. 


| C ¢ 
c 


As a matter of notation we write 
ri(we +k)? = f(w) + f(-w) 
where 
f(w) = | dy/(w — iy)(y? — k*)'? = (w’ + k’)'” log (11 + O)/(1 — Bb] (5.2 
and 
t = (w+ tk)'?(w — tk)”. (5.3) 


Upon using the definitions of the arg(w + ik) and arg(w — 7k) which we wrote above, 
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we find that the cut w plane (Figure 2) maps into the upper half of the ¢ plane by means 
of (5.3). Further the transformation z = (1 + #)/(1 — ®#) is bilinear and hence maps the 
upper half plane Jmt > 0 into the upper half plane Imz > 0 (Figure 2). We therefore 
define the branch of the logarithm by 0 < arg z < x. Given this information, we may 
verify that f(w) is regular in the lower half plane v < k. The factor (w? + k*)~'” multi- 
plying the logarithm in f(w) changes sign as we cross the lower branch cut. But the 
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Fie. 2. 
arg z = 0 on this cut so that the real part of f(w) is continuous. In crossing the branch 
cut log | z | becomes — log | z | so that the imaginary part of f(w) is continuous. Hence 


{(w) is regular in the required lower half plane. On the upper branch out the arg z = 7 
and hence the sign of the radical multiplier introduces a discontinuity. 

Similar comments apply to the function {(—w). The transformation t’ = (w —‘ik)'?- 
(w + ik)~'” sends the cut plane into the half plane Imt’ < 0 and the transformation 
2’ = (1 — tv’) (1 +”) takes the half plane Jmt’ < 0 into the upper half plane Jmz’ > 0 


with the upper branch cut mapped into the positive z’ axis. We put 0 < arg z’ < x to 
define the branch of the logarithm in f(—w), that is 
f(—w) = (w + kb)? log 1 — 011 + 


and we find that f(—w) is now regular in the upper half plane v > — k. 
As a check, we have in the strip — k <v <k 


f(w) + f(—w) = (w’ + k*)""”? [log z + log z’ 


= (w’ + k’)'”* [log | zz’ | + i argz + i argz’]. 
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Since z = —1/z’, this last expression reduces to the one with which we started. Hence if 

we write exp [a(w’ + k’)'’’] in product form as we have already described, we find that 
exp [a(u? + k’)'”?] = N_/N, 

where 


N_(w) = exp [a(w’ + k’)f(w)/z7] 


1/N.(w) = exp [a(w* + k°)f(—w)/ri]. 


The asymptotic forms of f(w) and f(—w) will be required at a later point in our 
work and since the analytical machinery is now developed, we shall provide these forms. 
Note first that 


log (1 + (1 — 2d) = log {i[w + (w* + k’)'”)/k}. 


“= < Oand — 2/2 < argw < 0. 


Further for u > 0,v < — k, — 2/2 < arg(w’ + k°) 
Hence for |w|—-o0,v< —k,u>O, 
log {7[w + (w? + k*)'”"]/k} = log [(2iw/k) + (ik/2w)] 
up to terms of 0(1/w). This in turn has the expansion 
(log 2iw/k) — k’/4w’. (5.4) 
Now the choice of the definition of the logarithm is dictated by the fact that the argument 
of the logarithm in (5.4) lies between the limits 0 and z. Hence if we consider the principal 


branch of the logarithm of w we see that log 7 = 7 r/2 and hence 
log (1 + O(1 — 7’ = in/2 + log 2/k + log w + O(1/w”) (5.4’) 


in this quadrant. If we now turn to the third quadrant, we see that we get certain minor 


< 2/2 


changes but the same dominant term. In this quadrant 0 < arg(w + k°)'* < 


while — 7 < argw < — 7/2. Then 
log (1 + (1 — t7' = log [ifw + (uw? + k*)'?}/k] 


= —log 2iw/k + O0(1/u") 


foru <0,v < —k,|w|—o. But again the logarithmic term has an argument between 
0 and x and since the arg w lies between — 2 and — 7/2, the negative of this argument 
lies between 7/2 and z. Hence as we might have surmised 


log (1 + )(1 — 27’ = —in/2 — log w — log 2/k + O(1/w’) 


for|w|—-o,v < + k,u < 0. This change in the asymptotic form of the logarithm of 
log(1 + #) (1 — 2)” in the third quadrant does not change the dominant term in the 
asymptotic form of f(w) however and (5.4’) is still in force. Similarly 


f(—w) = —|log w + log 2/k — ix/2 + O(1/w’)]/w 
for | w|—-«,v > — k. From this we conclude that 
N_(w) = exp [aw/2 + (aw/zi)(log 2/k) + (aw/mt) log w] 
if we retain the non-vanishing terms for | w | ~«,v < + k, while 
1/N.(w) = exp [(—aw/2) — (aw/xi)(log 2/k) — (aw/mi)(log w)] 


for|w|—-o,v> —k. 
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The factoring of 1 — 8/(w* + k’*)'” into the product of two terms, one of which is 
regular in the lower half plane v < 0, while the other of which is regular in the upper 
half plane v > — k is guaranteed by the original work of Wiener and Hopf. However, 
rather than apply their original methods, we shall follow an approach which will enable 
us to express the factors in terms of f(w) and f(—w). Again following the lead from 
the form (b), we write 


1 — 8/(w* + k*)'” = exp [In {1 — B/(w* + KY} ] 


where we choose the branch of (w* + k’)'”? which reduces to k for w = 0 as well as the 
principal branch of the logarithm. Instead of decomposing the logarithmic factor ad- 
ditively, we shall decompose its derivative in this fashion since differentiation will not 
alter the open strip of regularity. Hence we shall examine first 


d 


% In [1 — B/(w? + k*)'?] = wB/(w? + k’)[fw? + kY'? — Bl 
ao 


are ' eee : oa omy ~—; — —— ars 172 * a 3 ea (5.5) 
2(w — ik) 2(w + 7k) w—o Bljw +k} Biw —o} 

Now of these five terms, we can state the following. The first and third are regular in 
the lower half plane v < 0, while the second is regular in the upper half plane v > — &k. 
The fourth and fifth terms are regular only in the strips — k <v <kand-—-k<v<0 
respectively. The fourth term however is the derivative of —(w* + k’)'”?/8 which has 
been discussed in (b). The fifth term may be decomposed as follows: 


Ca so 
Biwe—o} 28 wto'w-—da)’ 

We know how to decompose (w® + k’)'” additively so that this last term breaks up 
into four terms, two of which are regular in the lower half plane v < 0, while the other 
two are still regular in the strip — k < v < 0. But the simple poles in these last terms 
may be extracted quite readily so that we finally have that the portion of (5.5) which 
is regular in the appropriate lower half plane is 

1 w 1 d 


ao? fe ee ae: Sah. eee ae © 2 2 k 
va 2(w — tk) si w—o rpidw w + K)f@) 


(w? + B)f(w) , (w* + B)f(w) 


2Bri(w + a) 2Bri(w — a) 


Ohta. , thi 


+ 2ri(w + a) 2ri(w — o) 


On the other hand, the term 


: 1 L @ oes ere 
U,(w) = - w + ik) - i du (w +k )f(w) 
I ee 
+ Mailw + 6) (we + k)f(—w) — Bf(+o)] 
] 


+e a [(w* + k*)f(—w) — B°f(-o)] 
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is regular in the appropriate upper half plane. Hence 


1 — B/p = | exp U_(w) aw || ex | U ,(w) aw | 


where the first factor is regular in the appropriate lower half plane while the second is 
regular in the appropriate upper half plane. The constants of integration are chosen so 
that the left side matches the right side in the strip of regularity. 

The asymptotic forms of exp {° U'_(w)dw and exp J” U,(w)dw are readily determined 


from our results in (b). We have for | w)>—-«0,v <0 
U_(w) = 1/w + 0(1/w") 
while 
Uw) = —1/w + O(1/u") 


for|' w|—-o,v> — k. Hence 
| exp U'_(w) dw = cw 
and 


exp ( (U) dw = C; ‘Ww 
in these respective cases where the constants ¢, and c satisfy the relation cc, = 1. The 
precise values of the constants do not concern us here. 

Now we may put together the factors (a), (b) and (c) to determine A_(w). In order 
to obtain an entire function of separation which is of algebraic growth, we introduce 
into K_(w) a factor exp x(w) where x(w) is to be chosen to make K_(w) of algebraic 
order for v < k, , | w| ©. This process has already been described elsewhere (1). We 
have then, forv < k,,|w{|—-« 

K (w) = 0(w"'”*) exp [(¢aw/r)(y — 1 + In ak/2r) + x(w)]. 


From this we see that if x(w) is chosen as 


x(w) = [Ll — y — In (ak/2n) |(taw/z) 
K_(w) = 0(w~’”’) in this lower half plane. Upon examining K,(w) for |w|—-2,v > —k 
we find that K.(w) = 0(w*'’*). Hence the factoring of K(w) and the asymptotic forms 


of the factors are known explicitly. 

In order to determine E(w), the entire function of separation, we merely apply an 
extended form of Liouville’s theorem. For example, for|w|—*, v < 0, E(w) = 
0(w~"), while for | w|o©,v > — k, E(w) = 0(w"””) as we can readily see from Eq. (4.4). 
From this we see that since E(w) is a polynomial, E(w) is identically zero. 


Earlier in our work we had assumed that ¢,(z, 0) was integrable from xz = 0 to any 
finite x. In view of the fact that J(w) = O(w™') and J(w) = O(w’’”’”) for | w | > @, v 
suitable limited, we see that J(x) is bounded for x — 0°, while ¢,(z, 0) = O(27'”*) for 


xz — 0°, thus verifying the integrability for finite z for these functions. The precise forms 
of ¢ will be investigated in the next section. 
VI. The Determination of ¢(x.y). In order to determine ¢(z,y) we return to the 
Helmholtz representations for ¢ in RP, and R, . For R, we have from (3.5) 
o(x, y) = [eosh any/ally Exp (2k2) + Y2 EXD | —ixx) | 
l 


—_ | GCG? (w, y, 0) L(w) exp (iwa) dw. (6.1) 
2r J, ; 
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Here G‘” (w, y, 0) is the bilateral Fourier transform of G"” (zx, y, 0, 0) an analytic func- 
tion of w in the strip — {k° + aj/a*}'” < v < 0. The convolution is justified here since 
L(w) is the transform of L(x) exp(vx) which is L(— ©, ~) for —-k <v < OandO <y 
< a (6). The path C is drawn in the strip of analyticity — k < v < O and the evaluation 
of the integral depends on whether x $ 0. If x < 0, we close below by a sequence of semi- 
circular ares passing between the poles of G” (w, y, 0) on the negative imaginary axis. 
In this case the integral becomes 

= [ Es Ky 2K. —- ale cosh p(y — a) + Bsinh ply — 2) lew (iwx) dw 

i(w — xk) i(w + x) K_(w)|p sinh pa — 8B cosh pa] p 


- Fees nK(— 2] exp (w,x)(o, + a°6") exp (w,) 


Qn J. 





<< |(—iw, — x) ° (x — iw,)_| iK_(—iw,)aw,(aB — a°B’ — a;) 
where w, = {a,/a° + k*}'”. The closing of the path C in the manner described above 
3/2 


is permitted since the integrand is 0 [w exp(— va)| for | w|—-o,v < 0. 

For x > 0, the path C in Eq. (6.1) is closed above. The singularities in the upper half 
plane are a branch point at w = 7k and four poles at w = + o and + mz. We introduce 
the positive imaginary axis from ik toi © as a branch cut. The phase of p has already 
been discussed on this cut, so the closing of C implies quarter circles on each side of the 
branch cut as well as integrals along the branch cut. This then supplies us with a path 
along which the integrand is single valued and analytic and Cauchy’s theorem may be 
applied for the evaluation of the integral. When we do so we get the residues at the poles 
+ o¢and + « and an integral along the branch cut. To do this we replace 

r exp (ap)(1 — B/p)K.(w) 
K_(w) = —*, $$ 
psinh pa — 8 Cosh pa 
and the integral along C becomes 
7 L i [ neko ‘ ralC.(— 2) || cosh p(y — a) + Bsinh Ay — «| ey 

2r Jo Li(w — x) i(w + x) (exp ap)(p — B)K.(w) 

The contributions from the quarter circles drop out as their radii becomes infinite and 
we are left with an integral along the cut as well as four residues. On the right side of the 
cut arg = 2/2 while on the left side of the cut arg p = —7*/2. Hence 


I + I = > residues 


or 
[ = —[cosh aoy/a][y, exp (ixx) + y2 exp (—ixz)] 
: Ress 4 2K (= 2 6° exp B(y — 2a) exp (ioz) 
o— kK o+k cK .(c) 
4 | 1K (0) mK 9 | B* exp B(y — 2a) exp (~iz) 
—(o + x) kK—-@o cK ,(e) 
~£ f° [nELo 4 nE co Tem Ty — ot pan Ty — a) 
xJ, L(w — x) (iv + x) K (iv) + o°’) 
k cos Ta + T sin Ta 
, azz. £19 FN 
exp (vz) 


where 7 = (v° — k*)'” 
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Now we examine ¢(z, y) in R, . Here we use the representation 


d(z,y = - | G(x, y, x’, 0)d,-(x’, 0) dx 


2 1 [ exp (twx) exp (py) L(w) dw 
p 


/ dw exp (wx + py) | x (x) 4 nk (0 | 
pK _(w) (w — k) (w + x) 


and again we distinguish the two cases, that is x 2 0. For x < 0, the only singularity in 
the lower half plane is the branch point at w = — ik. Hence we introduce a branch cut 
along the imaginary axis from — ik to — 7 © and choose a closing path of the variety 
we have just described for 0 < y < a, x > 0. Upon doing this we find that for y < 0. 
z<0 
d(z, y) = . | aoe ear | 11K +) + a Ca] cos Tydv,. 
rT. TK (—ww) | (—w — x) (—w + x) 

For y < 0, xz > 0, we close above in a fashion similar to the case x < 0, y < 0 to obtain 
a function ¢(z, y) which has the same functional form as we did for the region x > 0, 
0 < y < a. This incidently establishes the continuity of ¢(z, y) for x > 0, y = 0: 

VII. Remarks. We have now formulae for the region x < 0,0 < y < a,forz > 0, and 
for x < 0, y < 0. For each of these expressions the integral or infinite sum, as the case 
may be, is uniformly and absolutely convergent. Due to the presence of the favoring 
exponential when x + 0, we may carry differentiation under the integral sign (or inside 
the summation). It is possible, therefore, to verify directly that the Eq. (2.1) and the 
boundary conditions (2.1’), (2.1) appropriate to the region are satisfied. For x < 0, 
0 < y < a, we note that the expansion coefficients in the series are 0(n-**) for n > @. 
In view of the fact that we have decaying exponentials in this sum, it is clear that the 
asymptotic requirements of the region are met. 

Consider the integral in the representation for ¢ when xz > 0. We write its integrand 
as [exp vz] [h(y, v)]. Now h(y, v) is non-singular in (— ©, — k) and O(v *”*) independent 
of y as v ~—-, It follows e** and h(y, v) are L*(— ~, — k) integrable and we have 


by Schwarz’s inequity 
— j ose 1/2 
| e*h(y,v) dvi < | | e”*” dy | h(y, v) |° dv 
. k e « k 


< al | r iv| — 0 asx — @ independently of y. 


This shows the asymptotic requirements far to the right of the barrier are met. 

The Riemann-Lebesque lemma shows the integral in the representations for ¢ when 
x < 0 or when z > O goes to zero as y becomes negatively infinite. Since this is true of 
the exponential terms in the case x > 0, we see that ¢ decays with depth. Clearly then, 
all of the order conditions we have had to assume in the formulation of the representa- 
tions of d(x, y) for R, and R, are satisfied. 

VIII. Reflection and Transmission Coeffcients. From the results which we derived in 
Sec. VI, we may supply the reffection and transmission coefficients of the traveling 
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surface waves for y = a,x 2 0. We note that for x -~@ 
d(x, ¥) = lys exp (tox) + ys exp (—tox)] exp (By) 
and this in turn is equal to 


-|? K .(k) WK (- | B° exp (8)(y — 2a) 
Go 


- exp (102% 
— kK o+k o K (a) P (toz) 


sa | vii (x) 4 YK é K) exp (—ioz). 


| 8° exp B(y — 2a) 
—ie + K a 


o K .(-—«o) 


From this we have then 





ae _[ nk.) mK (— 9) 8° exp (—28a) (8.1) 

ro lo—k o+k o K (oe) maa 
and 

ws ly K(x) yoK C9] 8° exp (— 28a) : 

ro ere ™ o—-k ¢ K.(-—o) ~ (8.2) 


That is, we have two linear relations between the amplitudes of the traveling waves in 
the regions x > 0 and x < 0, y = a. The conditions (8.1) and (8.2) are independent as 
one can readily see by evaluating the determinant of the system (8.1) and (8.2). 

As usual we have a relation between the magnitudes of y, , y. , y3 and y, . We form 
the integral 


| (0 - 0%.) as! (8.3) 


which is equal to zero since ¢ is source free and sufficiently regular. This integral is 
evaluated over the “boundary” of the channel with due recognition for the barrier. In 
view of the particular asymptotic forms we have found for ¢ for | z | ~#, v ~—© 
(8.3) reduces to 


v2 —lyw P= En? — lve |?) 


Bax(Ba — B’a? + ao)[exp (— 28a) ][1 + exp (—2a,)] 
SS ee eee, (8.4) 
20a,(a@) — ap) 
The expressions we have found for 7; and +, satisfy this relation. On the basis of (8.4) 
we can define the reflection coefficients and traismission coefficients. Note that if y, = 0, 
we have transmission to the right and reflection to the left. | y. | = | y2/y, | is the magni- 
tude of the left reflection coefficient. Then 


¥3 |" Qeai(ao — af) exp (2a) 





vy, | Bax(Ba — B’a" + a)(1 + exp (—2a)] 


defines the transmission coefficient to the right, that is | ¢e |. We find in this case 
= & 9 
Int = (S54) 
17. | a + K ’ 


| tr |? = L .) 


(¢ + x)” 
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Similarly | rz | = | rz! and | t; | = | te | , that is for the case in which y, = 0. Note that 
for the case k = 0 
; —_— R oi (28 —22)) 
sind Ene i Makan aB + a/’ 
jw 14 Po 4aBao ' 
(aB + a) 


Had we only been interested in magnitudes of reflection and transmission coefficients, 
these results might have been anticipated from classical transmission line theory. How- 
ever, the determination of the various reflection and transmission coefficients involves 
a phase for each of these coefficients. These angles may be calculated from K,(7), 
K. (—x), K. (c) and K, (—c) but we shall not pursue this point here. 
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—NOTES - 


THE FUNDAMENTAL THEOREM OF ELECTRICAL NETWORKS* 


By J. L. SYNGE (School of Theoretical Physics, Dublin Institute for Advanced Studies) 


I wish to thank Messrs. L. C. Robbins, Jr. and W. K. Saunders for drawing my 
attention to what may appear to be a flaw in the reasoning of the paper with the above 
title." The proof there given is not convincing, and I take this opportunity to supply a 
fuller proof, in which the reasoning is a little delicate, but which seems to establish the 
result. 

Consider the two propositions: 


(A) e=0 implies i= 0. 
(B) Z’ ~ 0. 


It is stated on p. 127 of the paper cited that (A) implies (B); that is the result we have 
to establish here. 
We have the following equations: 


fer 


(1) i = Ci’, e’ = Cre, e’ = Z’i’. 


These equations embody all our knowledge about the behavior of the network, and any 
set of vectors i, i’, e, e’ consistent with them are to be regarded as possible. But we 
must bear in mind the definition of the mesh currents i’ as branch currents (components 
of i in branches-out-of-tree); this means that the matrix C is such that 


(2) i’ ~ 0 impliesi + 0. 
Choose e = 0 and e’ = 0. Then (1) read 
(3) i = Ci’, 0 = 0, 0 = Zi’. 


Any vectors i, i’ satisfying these equations are possible. 

Suppose that (B) is false, i.e. suppose Z’ = 0. Then there exists a non-zero i’ satisfying 
the last of (3), and, by (2), the i given by the first of (3) is non-zero. Thus, on the as- 
sumption that (B) is false, we have a solution of (1) with e = 0 and i ¥ 0. Hence, if 
(B) is false, (A) is false. Therefore if (A) is true, then (B) is true. In other words, (A) 
implies (B), which is the required result. 


*Received June 1, 1952. 
1Q. Appl. Math. 9, 113-127 (1951). 
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PARALLEL REFLECTION OF LIGHT BY PLANE MIRRORS* 


By JOSEPH B. KELLER (Jnstitute for Mathematics and Mechanics, New York University) 


1. Introduction. If three plane mirrors are mutually orthogonal then, as is well 
known in geometrical optics, any ray’ incident on these mirrors in the octant which 
they bound will be reflected once from each mirror and will emerge parallel to its original 
direction. Similarly any ray incident on two orthogonal plane mirrors, in the quadrant 
which they bound and in the plane perpendicular to them, will be reflected once from 
each mirror and will emerge parallel to its original direction. It is interesting to investi- 
gate whether there are other configurations of plane mirrors which also have this property 
of parallel reflection. 

For two mirrors, with the incident ray in the plane perpendicular to them, the answer 
is easily found although it does not seem to be well known. It is this: if and only if the 
mirrors meet at an angle 7/n, where 7 is any even integer, every ray will emerge parallel 
to its original direction. Furthermore every ray suffers n reflections; for example, two in 
the case of orthogonal mirrors discussed above. 

For three mirrors meeting at a point Synge” has shown that if every ray is reflected 
exactly once from each mirror every ray will emerge parallel to its original direction if 
and only if the mirrors are mutually orthogonal. However, in view of our above result 
for two mirrors, it is reasonable to drop the restriction that every ray be reflected exactly 
once from each mirror. We then find that there are indeed other configurations of three 
mirrors meeting at a point which yield parallel reflection for every incident ray. These are 
the configurations in which the sets of dihedral angles are (4/2, 4/3, 2/4), (4/2, 2/3, 2/5) 
and (7/2, 7/2, r/n) where n is an even integer. The numbers of reflections suffered by 
any ray are 9, 15 and n + 1 respectively. The case (1/2, x/2, r/n) with n = 2 is the case 
of mutually orthogonal mirrors, in which case each ray is reflected three times. Further- 
more, these are the only configurations of three or more mirrors meeting at a point which 
have the parallel reflection property. 

In section II we formulate the problem, and in section III we deduce the result 
described above. 

2. Formulation. Suppose three planes meet at a single point and that one side of 
each plane is designated as the reflecting side. Then one and only one of the eight regions 
into which the planes divide space will be on the reflecting side of all three mirrors. Any 
oriented straight line in this region which can be extended indefinitely in its negative 
direction without intersecting a plane is called an incident ray’. If an incident ray inter- 
sects a plane it is reflected according to the law of reflection; the reflected ray may be 
again reflected, etc. If a reflected ray can be extended indefinitely in its positive direction 


*Received June 2, 1952. This work was performed at Washington Square College of Arts and Science, 
New York University and was supported in part by Contract No. AF-19(122)-42 with the U.S. Air Force 
through sponsorship of Geophysics Research Division, Air Force Cambridge Research Center, Air Ma- 
teriel Command. 

1Rays which intersect an edge or vertex or are paralle] to one or two mirrors are excluded. The latter 
rays do emerge parallel to their original directions although they suffer only two or one reflections, re- 
spectively, not being reflected from the mirrors to which they are parallel. 

2J. L. Synge, Reflection in a Corner Formed by Three Plane Mirrors, Quarterly of Applied Math., 4, 


166-176 (1946). 
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without intersecting a plane, it is called a final reflected ray. The problem we consider 
is that of finding all configurations of three planes meeting at a point with the property 
that every final reflected ray is parallel to the incident ray from which it originates. 

First, following Synge, we shall reformulate the problem in terms of reflection of 
points on the surface of a sphere. To this end we construct a sphere with center at the 
meeting point of the planes. The surface of the sphere intersects the planes in great 
circles. The trihedral region on the reflecting sides of all three planes cuts out on the 
sphere a spherical triangle 7’, which is bounded by the three great circles. We represent 
a ray by the point of intersection, with the sphere, of a parallel infinite line through the 
center. Since a straight line through the center intersects the sphere twice the more 
negative of the two points is chosen. Thus all parallel rays with the same orientation are 
represented by the same point, while antiparallel rays are represented by antipodal 
points. It is clear that all incident rays, defined above, are represented by points in the 
spherical triangle 7’. 

Now suppose that a ray, represented by the point P, is reflected in a plane, represented 
by the great circle A. Then the reflected ray is represented by a point P, which is obtained 
by reflecting P in A. This can easily be seen if the point of reflection is transported to 
the center of the sphere. The ray from P to the center is the incident ray; the ray from 
P, to the center is the negative extension of the reflected ray. 

Since, by definition, only one side of each plane is reflecting a point P must lie on 
the reflecting side of a great circle A in order to be reflectable in A. We will call a re- 
flection of P in A admissible if P does lie on the reflecting side of A. The reflecting side 
of A is the side or hemisphere containing the spherical triangle 7’, according to the 
definition of 7. Now a ray emerges, or is a final reflected ray in accordance with the 
definition above, only when it is not reflectable in any plane. Thus a point P represents 
a final reflected ray only if it lies on the non-reflecting side of each plane. But every 
point which lies in all three hemispheres not containing 7 must lie in the spherical 
Triangle 7’ antipodal to 7. We can now formulate our problem as that of finding all 
spherical triangles 7 with the property that every sequence of admissible reflections of 
any point P in T leads to the antipodal point P’ in T’, since P’ represents a final ray 
parallel to the incident ray represented by P. 

3. Solution. Let us suppose that 7 is a triangle with the desired property and that 
P, and P, are two points in T. We first want to prove that if an admissible sequence of 
reflections R, carries P,; into some point P, and an admissible sequence of reflections 
R, carries P, into the same point P, then P, and P, are identical. In other words, no 
two distinct points can have a common image under admissible sequences of reflections. 
To prove this assertion we note that, by hypothesis, every sequence of admissible re- 
flections of a point of 7 terminates on the antipodal point. Thus, in particular, the 
sequence R, which carries P, into P can be continued by a sequence of admissible re- 
flections which we shall call S, so that R, + S, carries P, into P{ , the antipodal point. 
But then S, carries P into P{ . Therefore R, + S, carries P, through P into P{ . Since 
this is an admissible sequence of reflections, it must by hypothesis terminate on P3 . 
Therefore P{ is identical with P; and so P, is identical with P, . 

Now let us consider the neighborhood of a vertex of T of vertex angle a. We wish to 
consider successive reflections of a point of T in the two sides meeting at the vertex. 
We introduce polar coordinates with the vertex as origin and one side as polar axis. If 
¢ is the polar angle the other side is given by ¢ = a. If ¢ = 6 is the angular coordinate 
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of a point of 7, its image after n successive reflections in the sides 0 and a will have the 


coordinate ¢,(@) given by: 


na + 6 n even | 
¢,(0) = if the first reflection is in side 0, | 
—(n — lha — 6 n odd 
(1) 
| 
(n — lha — 6 n odd 
¢,(6) = if the first reflection is in side a. | 
n even 
—na + 6 J 


These reflections are admissible only until the image falls in the vertical angle x to r + 
a (or — 7 to — + a) since then the image is on the non-reflecting side of both mirrors 
meeting at the vertex. [To prove the result stated previously for two mirrors, we set 
¢, = = + Gor — r+ Gand find a = 2/2 with n even.| 

From the expressions in Eq. 1 it is not difficult to show that some image of every 
point will lie in the vertical angle 7 to + a. Let us consider the least value of n for 
which the image of any point of the original angle lies in the vertical angle. The sequence 
of admissible reflections which carries a particular point into the vertical angle in this 
number of reflections is necessarily an admissible sequence for all points in the original 
angle. Therefore the image of the original angular region under this sequence is another 
angular region which overlaps the vertical angle. If this overlapping is complete, then 
one may deduce, as above, that a = 2/n where n is an odd or even integer. If the over- 
lapping is incomplete the image covers an edge of the vertical angle. The points outside 
this edge may then be reflected in it, and some of them will fall on already covered points. 
Therefore distinct points of the original angle will have a common image, contradicting 
the theorem above. Thus the overlapping must be complete and therefore, we conclude 
that each vertical angle is of the form 2/n, where n is an integer, not necessarily the same 
for all vertices. 

Combining this result with the fact that the sum of the interior angles of a spherical 


triangle exceed 7, we obtain the inequality. 


T T T 
74747510. (2) 


Ny Ne Ig 


The only solutions of this inequality are: 
\ 


(z m =) (z x *) (z =) (eee (3) 
oe ig 2’ 3’ 4/’ 2’ 3’ 5/’ 7 na)” 


Thus the only possible triangles having the desired property are included here since 
these solutions were obtained by imposing necessary conditions. 

By examining the admissible reflections of these triangles on a sphere, we find that 
only the following actually do yield parallel reflection: 


(z£2) z,¢ 2) Es ¢) = ' 
9’ 3’ 4)’ 2’ 3’ 5 . 22’ n ~n even. (4) 


The other triangles do map onto their antipodal triangles but with the sides interchanged; 
in other words, a point does not land on its antipode but on the antipode of its image in 
the principal altitude of the triangle. In the same process we find that the respective 
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numbers of reflections are 9, 15 and n + 1. This completes the proof for the case of three 
mirrors. 

For more than three mirrors all the considerations preceding equation 2 are valid, 
but the corresponding inequality has no integer solutions at all. 

To prove only Synge’s result, which is included above, a much simpler method 
suffices. If a point arbitrarily near a vertex is reflected successively once in each of the 
sides meeting there the image remains arbitrarily close to the vertex. The third reflection 
must bring it to the antipode, an angular distance arbitrarily close to 7 radius. Thus 
the angular distance from the vertex to the third side is arbitrarily near 7/2, and this 
applies to each side. Thus the triangle must be an octant of a sphere, which completes 
the necessity proof. The sufficiency is proved by reflecting an octant three times in its 
sides. 

It is of interest to compare our results with those described on pps. 81-83 of Coxeter’s 
“Regular Polytopes’”’ (Methuen, London, 1948). There the discrete groups of reflections 
generated by three or more concurrent planes are deduced. It is first found that all the 
dihedral angles must be of the form 2/n, then equation 2 is applied, and the solutions 
are given by equation 3, with none for more than three planes. In the problem of discrete 
groups all sequences of reflections are permitted, while in the present problem only 
‘‘admissible’”’ sequences of reflections are permitted, which is a weaker condition. On 
the other hand we have the condition that some admissible sequence must carry each 
point into its antipode, and this condition is strong enough to limit the solutions to those 
given in equation 4—a subset of the configurations. generating discrete groups. 

The new parallel reflecting configurations given herein may have practical value in 
some applications because they reflect back only those rays which are incident from a 
certain angular region, and the angular region may be made much smaller than that of 





the mutually orthogonal mirrors. 
I wish to express my indebtedness to my colleague, Prof. Wilhelm Magnus, who 


supplied a crucial part of the necessity proof giv 7e. 
pplied l t of th ty proof given above 


NOTE ON SELF-PROPAGATION OF TURBULENT SPOTS* 
By CARL E. PEARSON (Harvard University) 


Introduction. Emmons (1), (2) has recently directed attention towards the manner 
in which transition occurs from a laminar to a turbulent boundary layer. It appears that 
the dominating phenomenon is the spontaneous generation of a large number of tur- 
bulent ‘spots’ which grow rapidly and eventually coalesce. It has been observed ex- 
perimentally that, except perhaps immediately after generation, each spot grows in 
such a manner as to maintain geometric similarity. Because the spot is growing in a 
boundary layer, the rate at which an elemental portion of its periphery propagates its 
turbulence outwards is dependent on the orientation of that portion; this dependency 
of propagation rate on orientation is necessarily associated with the shape of the spot. 

It is the purpose of this note to prove the simple but interesting result that for either 
of two reasonable speculative hypotheses concerning propagation rate, the shape of the 
spot approaches a certain asymptotic shape which is independent of the initial spot 
shape. 


*Received July 18, 1952. 
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Material point hypothesis. It is possible to conceive of the propagation of turbulence 
into a non-turbulent region as being due to the effect of small random eddies being con- 
tinually induced in the quiescent fluid immediately adjoining the turbulent region. In 
a laminar boundary layer, the propagation of such induced turbulence may be con- 
sidered as being at a rate which is dependent on the orientation of the turbulent surface 
in question and in a direction perpendicular to the surface; this leads us to the following 
equivalent picture. 

Suppose that at every instant of time a set of material points occupy a closed convex 
curve in the plane. Assume that each material point moves perpendicularly to the curve 
on which it lies and with a velocity R(6@) depending only on the slope @ of the curve at that 
point. Let the family of curves generated in this way be represented parametrically in 
the complex plane by 
2(6, t) = x(8, t) + 2y(8, 0) 


i being the time. Consider first the rate at which the @ value associated with a particular 


material point changes in time. 
In Fig. (1), let A and B be positions occupied by the curve at times ¢ and ¢ + Al, 











Fia. 1. Neighboring trajectories. 
and let C and D be two neighboring orthogonal trajectories (i.e., paths followed by two 


neighboring particles). Then to the first order, the change in 6 for a particle on trajectory 
D occurring in a time interval At is 


- (#40 31) 
Aé = (fi 66 At 0s, 


dR 66 = R(6 + 66) — R(8), 


dé 
a0 _ ak fds 
dt dé/ dé 


Now the velocity of the material point is 


where 


and therefore 


dz 0z dz dé 


di at* aode 
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whence 
—_—. (ahide = )(22 2) “5 (ae ) 
at = — te” + Var 7a6 ° "Nag + 30) = © \ao — **)- 
But the right hand side is a function of 6 alone, and the equation may therefore be 
integrated to yield 


(dR. 

z= ot _— iR)t + C(6), 
dé 

where C(@) represents the parametric position of the curve at time ¢ = 0. Thus the curve 


approaches the asymptotic motion 


r= nt Rsin 6+ a cos ), 


2 
y= (a! sin 6 — R cos 0). 
(These two equations may also be used to determine R is the asymptotic shape is known 
from experiment). 

Wavelet hypothesis. Each point on the periphery of a turbulent spot may be re- 
garded as the continual origin of a disturbance which spreads outwards from it in time; 
the instantaneous envelope of these disturbances is then the self-propagating curve 
itself. The shape of the disturbance surrounding any point of origin need not be circular 
if the region into which it spreads is inhomogeneous (as in the case of a laminar boundary 
layer). Let the velocity of the (infinitesimal) disturbance be R (8) as shown in Fig. (2). 








Fic. 2. Wavelet velocity. 


Then if the equation of the closed curve is 
x = f(6, 2), y = g(, 0), 
the equation of a wavelet emanating from the point (6) will be, after time At, 


zr— f(0, t) AtR(8) sin B, 


y — g(6, t) = AtR(8) cos B. 
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The equation of a neighboring wavelet will be 
xz — f(0@+ Aé, t) = AtR(B + AB) sin (8 + As), 


y — g(0+ Ad, d) — AtR(B + AB) cos (8 + ABs), 


Il 


and hence the relation between 6 and £ is obtained from 


a Ad = —Al < (R sin B) AB, 
og A@ = At : (R cos B) Ag, 
06 0B 
whence 
09/00 _ _ A(R cos 8) /d8 


tan 6 = : ; aoe (1) 
' O(R sin 8)/dB 


Of /d0 in 
This equation gives the direction 6 of the intersection point of the wavelets as a function 
of @. 
Now consider again the curve to be composed of material points, this time moving 
in the §-direction of Eq. (1). 


1 Jy(6 + Ad, t+ At) — y(O, t+ Al y(0 + A@, t) — y(8, bt) 


d 
— (tan @) = l 
at ‘ wake ao s1-0 At \x(6 + AG, t+ At) — 2(0,t + Ad x(6+ dB, t) — x06, 0) 


li ] fay 060 — At(dB/dé) A(R cos B)/dB Oy 0\ 
= lim ) — oe 
aro At \dx/00 + At(dB/dé@) d(F sin B)/dB 0x /086) 
dy/06 f(de dé) A(R cos B)/dB ‘ (dB/d8) a(R sin B)/dB 
62/88 | 04/06 0x/06 
dg/de f a ,) ais 
== {. (R cos B) + tan 8 : (R sin 8B)? = 0, 
02/06 \dB op ) 
by (1). Then 
dx ‘ Ox dx dB Or 
= Rsn 6 = — xe =; 
dt ot + 06 dt ot 
whence 
x = (Rsin B)t + F(6). 
and similarly 
y= —(R cos B)t + G(8), 
where F(8) and G(8) are the parametric representations of the original curve, 8 being 
now used as the parameter for convenience. Thus there is an asymptotic shape, in fact, 
that of Fig. (2). 


REFERENCES 
(1) H. W. Emmons, The laminar-turbulent transition in a boundary layer, Part I, J. Aero. Sci. 18 490 


(1951). 
(2) H. W. Emmons, Part II, to appear in the Proc. of Ist U. S. National Congress. Theor. Appl. Mech. 

















1953 M. G. SALVADORI AND F. DIMAGGIO 223 


ON THE DEVELOPMENT OF PLASTIC HINGES IN RIGID-PLASTIC BEAMS* 
By M. G. SALVADORI anpb F. DIMAGGIO (Columbia University) 


Synopsis. A study is made of the development of plastic hinges in free-free, rigid- 
plastic beams acted upon by certain types of distributed, dynamic loads, characterjzed 
by a concentration parameter c, whose value varies between zero and infinity as the load 
distribution varies from a uniform to a concentrated load. 

It is found that as the load intensity increases a plastic hinge is first developed at 
the center of the beam, whatever the value of the concentration parameter. 

As the load intensity increases beyond this value, the development of successive 
hinges depends upon the parameter c: two additional hinges develop laterally if ¢ > 3.45, 
while the central hinge first splits if ¢ < 3.45 before lateral hinges can develop. In the 
first case the lateral hinges wander toward the center of the beam. In the second a wider 
central portion of the beam becomes plastic. 

1. Introduction. The elasto-plastic behavior of free-free beams under dynamic 
loads or initial velocities has recently been studied by Bleich and Salvadori’, who 
obtained expressions for the elasto-plastic displacements and the plastic deformations 
in terms of infinite series expansions of normal modes. 

Bleich and Salvadori showed that upper bounds for the plastic angle developed at 
the center of a free-free beam under symmetrical loads or velocities may be obtained 
by considering the beam as a rigid-plastic body (as suggested by Prager for concentrated 
dynamic loads). These bounds are good approximations of the true plastic angle when- 
ever the dynamic loads or initial velocities are so high as to produce bending moments 
several times larger than the capacity moment. Another approximation of the plastic 
angle may be obtained by freezing the beam into a rigid-plastic body at the time when 
the first hinge is developed. In the case of initial velocities this approximation is often 
a lower bound of the plastic angle. 

Lee and Symonds’ obtained upper bounds for the plastic angle and studied the 
successive development of plastic hinges in rigid-plastic, free-free beams acted upon 
by a normal concentrated dynamic load P at the middle point. They proved that, after 
the appearance of a central plastic hinge, two additional hinges develop laterally, which, 
as the load increases, wander toward the center of the beam. 

In what follows the development of plastic hinges due to a symmetrical distributed 
load p is studied by means of a concentration parameter c, which permits the distribution 
of p to vary from a uniform load po, corresponding to c = 0, to a concentrated load P, 
corresponding toc = ©. 

2. Rigid-plastic beams. A rigid-plastic beam is an idealized beam of infinite rigidity 
capable of becoming plastic suddenly at a section where the bending moment due to the 
loads and the inertia forces reaches a given value, the so-called capacity value M, . The 
material of a rigid-plastic beam has therefore a stress-strain diagram of the type indi- 


*Received June 10, 1952. This paper presents results of research conducted at Columbia University 
under contract Nonr-266 (08) with the Office of Naval Research. 

1H. H. Bleich and M. G. Salvadori, Vibration analysis of elasto-plastic structures, with application 
to impulsive motion, Office of Naval Research Project NR-360-002, Contract Nonr-266(08), Technical 
Report No. 1, (1952). 

2. H. Lee and P. S. Symonds, Large plastic deformations of beams under transverse impact, J. Appl. 
Mech. 19, 308-314 (1952). 
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cated in Fig. 1. Let us consider, in particular, a rigid-plastic, free-free beam of length 
21 and constant mass m per unit of length, acted upon by a distributed normal load 
p(x, t) and referred to an z-axis with origin at its middle point (Fig. 2). Under the as- 
sumption that the load p be symmetrical with respect to the origin and maximum at 
the origin, the bending moment will at first be maximum at 0, since the average pressure 
P/(2l) and the constant inertia forces per unit of length are equal and opposite, and 
the moment of the pressure p(z) —P/(2l) is maximum at 0 (Fig. 3). In this first phase 
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of the motion the beam translates rigidly with a constant acceleration a in the y direc- 
tion and, as p increases in magnitude, say linearly, there comes a moment when the 
bending moment M at 0 reaches its capacity value M, . At this point a hinge is suddenly 
developed at 0 and each half-beam acquires an angular acceleration —a on top of the 
linear acceleration a. From this point on the moment of the rotational inertia forces 
contributes to the bending moment, which may reach its capacity value at sections 
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other than x = 0. An analysis of this second phase of the motion, limited to the develop- 
ment of additional lateral hinges or to the splitting of the central hinge, is contained in 


what follows. 
3. The concentration parameter. The type of load considered is distributed sym- 
metrically about the origin according to the exponential law: 


c 


lA 


pz) = pe a (1 + cele” 0s2zs)) (1) 


bo 


where: 


z= 2/l (2) 


and c is a constant called the concentration parameter. 
The total load on the beam due to the distribution (1) equals 


al 
P= pls 2 | (1 + cze “* dz 





(3) 
- 2p.i(1 mgt os ce) 
and the pressure p, is so chosen to make P independent of c, so that 
P/(2) 
p. = ic ° (4) 


1—e (1 + c/2) 
It is seen from Eqs. (1) and (4) that as c approaches zero p(z) approaches a uniform 
distribution p,) = P/(2l), while as c approaches infinity p(z) becomes a concentrated 
load P applied at « = 0. Figure 4 gives the graphs of p(z) for c = 0, 2, 6 and 10. 
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4. The rigid translation phase. The acceleration a of the beam in the y-direction, 
before the bending moment reaches its capacity value at 0, is given by: 


P 


i (p./m){1 — e “(1 + c/2)]. (5) 


a= 
The corresponding inertia forces —ma and the load p(z), per unit length, create a moment 
M(0) at the origin evaluated by the rotational equilibrium of a half-beam about 0 


(Fig. 5): 
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a P(2) 


| on fp! 7 P . v 
M(0) = —p * | (1 + cz)z *dz+ Pp D [ aa (1 +8] 


/? 
= Pe [1 —3/e+e (2+c¢/2+3 c)] (6) 


It is convenient to introduce a non-dimensional load parameter u defined by the equation 


__ es 7) 

es M, ’ ‘ 

where P is given by Eq. (3) and M, is the plastic capacity moment of the beam. In 
terms of u the moment 1/(0) becomes: 





=~ Sle + 0°82 + 0/2 +84 . 
— 3/e +e (2 + 6/2 + 3/e) | (8) 


M(0) = M op L : 
4fl — e-(1 +¢/2)] 


A central hinge develops and the rigid translation phase of the motion ends when M (0) 
reaches the value M, , that is for the value up of » defined by: 





{ 1—e (1 + c/2) oe (9) 
1—3/e+e°(2+c¢/2 + 3/0) 


It is seen from Eq. (9) that the value yo of u at which the central hinge appears approaches 
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4 as c approaches infinity (concentrated load). As c approaches zero yo becomes inde- 
terminate, but expansion of e* into a power series gives the asymptotic value 

weS-8 «1 (10) 
which approaches infinity as c approaches zero. The graph of uo versus ¢ appears in 
Fig. 6a. 
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Fic. 6a. 


5. The translational-rotational phase. For values of u larger than yo the half-beams 
translate with an acceleration a and rotate rigidly with an angular acceleration —a, so 
that the acceleration of the section x = lz of the beam becomes: 


a(z) = a — alz. (11 
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The two unknown accelerations a and al are determined by applying the equations of 
dynamics in translation and rotation to the half-beam (Fig. 7): 
ol 
P/2— | ma(z) dz = 0 (12a) 
ol al 


M,+ | p(zjz dz — m | a(z)z dz = 0. (12b) 


0 “0 
By means of Eqs. (11), (1) and (3) these equations reduce to: 
ma — $mal = p,{1 — e “(1 + c¢/2)] 


$ma — mal = M,/l° + 3p.[3/e — e (3 + 3/e + 0)] 
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and give the following values of ma and mal: 
ma —6M,/U? + p.[4 — 9/e +e°(5 + 9/e + 0)) (13a) 
mal —12M,/l* + p.|(6 — 18/ce +e “(12 + 18/e + 3c) | (13b) 
y 











plz) 


Fia. 7. 


By means of Eqs. (1), (7), (11) and (13) the bending moment at any section x = lz 
becomes (calling u the dummy variable of integration) 


M(iz) = | [p(u) — ma(u)|(u — z) du 


1 


y Mo fon — e“(1 +.¢/2)\(1 — 32° + 22’) 


2i1 — e “(1 + ¢/2)] | 
4 ; {(3/e — 2z) —e “(2 + 3/c) 
+ 2°[(4 — 9/c) +e°(5 + 9/e + 0] 
—2(2-—6/e+e (4+ ve + on}. (14) 
A second hinge will develop in the beam as soon as M(z) has a stationary value at, 
say 2 = z, * 0, and this stationary value equals the capacity moment M, in absolute 


value. Since M(z) is certainly maximum at z = 0, even after the hinge has developed 
at. 0, because of symmetry, the stationary value of M at z, must be a minimum (Fig. 8) 
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or M(z) would have a maximum larger than M, at another section. Hence the section 
z, , at which a second hinge develops as yu increases beyond yp to a value yp, , is defined by 
the conditions 

M(2) aM 

(15) 


Mt ae 


which by means of Eq. (14) become: 
2(1 — e “(1 + ¢/2)\(2 — 32° + 22") + w/2{(3/e — 22) — "(2 + 3/c) + 27[(4 — 9/0) 
+¢‘'(5+9/e+ 0] —2 [2 —-6/e+e€'(4+ 6/e+ 0]} = 0; (16a) 
12/1 —e “(1 + ¢/2)\(@ —2z) + w{-—[l — e (1 + 3cz)] + 2[(4 — 9/0) 
+¢(5+9/e+ 0] — 32° [(1 — 3/c) + e°°(2 + 3/e + c/2)]} = 0. (16b) 


Solution of the simultaneous Eqs. (16a) and (16b) for various values of ¢ gives the 
graphs of uw, and 2, versus c plotted in Figs. 6a and 6b, respectively. It may be shown 
analytically that as ¢ approaches infinity », approaches 22.9 and z, approaches 0.4; 
while, as c decreases, the second hinge moves outward and y, increases, approaching 
infinity as c approaches a value between 2.0 and 2.1. 

6. The splitting of the central hinge. The values of z, and y, obtained above for 
the development of a second hinge in the half-beam were predicated upon the fact that 
M(z) would remain a maximum at z = 0. This is guaranteed by the equation 


Fal! 
ER < 0. 


If instead, as uw increases above yo , the second derivative of M becomes positive for 
values of u such that 
Mo <b < M1, (17) 


the moment M would become minimum at 0. Since this cannot happen lest M(z) be 
greater than M, at another section, the central hinge splits into two symmetrical hinges 
which wander outward plasticizing a central portion of the beam, in which the moment 
has the constant value M, . The value uo,, of u at which the central hinge splits is there- 
fore characterized by the condition 


a°M 
—., = ( 
| Oz |. », (18) 
that is, by the equation: 


[4—9/e-—c/2+¢ (5+ 9/e+o]u + 6-2 + € (2+ 0)] = 0, 


from which: 





- seen ce e€ 2+o] “a (19) 
(-9+ 4c -—c/2) +e ‘(9 + de + Cc) 

The graph of uo,, versus c, plotted in Fig. 6a, has a minimum for c = 2.2, approaches 
infinity as c approaches zero and 3.63, and crosses the graph of uw, force = 3.45. 

It is thus proved that lateral hinges will develop in the beam when yu = yu, while the 
central hinge remains stationary if c > 3.45. If instead c < 3.45, the central hinge will 
split, plasticizing a wider central portion of the beam, when uw = wo,, , before lateral 


hinges may appear. 
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THE EFFECT OF COALESCENCE IN CERTAIN COLLISION PROCESSES* 
By Z. ALEXANDER MELZAK (McGill University) 


This note is concerned with the variable size-distributions of atmospheric and colloidal 
particles. The following physical model is adopted: a large number of particles are 
randomly distributed in space; these particles are in random motion due to e.g. Brownian 
effect or turbulence; the system of particles is given by the particle-mass distribution. 
The distribution varies as a result of collision processes in which the randomly moving 
particles meet and coalesce. The particles do not break up and do not enter or leave the 
system under consideration. It is assumed that the total number of particles is large 
enough to justify the use of 

(a) the mass density function f(z, 0); f(a, dx being the average number per unit 

volume of particles of mass (or volume) z to x + dz at time ¢, 

(b) the coalescence function ¢(u, v), where 

N(u, v, 2) = flu, Of, Deb(u, v) du do dt, (1) 
Niu, v, t) being the average number, per unit volume, of collisions resulting 
in coalescence, between particles of mass u to u + du and v tov + dv respectively, 

in the time interval ¢ to ¢ + dd. 
The function (u, v) takes care of such factors as: the collision cross-section, the mobility 
of the particles, the fraction of collisions resulting in coalescence and any shape factor 
affecting the collision. With these definitions, the condition for the conservation of mass 


becomes: 
Nu, v, ) — N(u, 2, 0). (2) 


“u=0 


— D dx dt = 5 
This means that the increase of the number of particles of mass x in time dt is equal to 
half the number of coalescences between particles of masses u, and z — u,0 < u < 2, 
minus the number of coalescences between the particles of mass z and any others. Sub- 
stituting equation (1) in equation (2) gives the fundamental equation of the process: 


“uteo= 


OO = S| iy, Of — y, O6y, 2 — dy - flz,) | fy, 00@ dy. ©) 
Cc 0 “0 
If N(é) is the total number of particles (per unit volume) then: 
N(t) = [ f(z, t) dz, (4) 
dN(t) _ ~ [ - . ¥ 


The initial condition will be taken as f(z, 0) = g(x), a known function corresponding 
to a known initial distribution. If this initial distribution is nearly homogeneous, and if 
¢(z, y) does not oscillate rapidly, then for small ¢ equation (5) may be approximated by: 


dN irr cN* _ 


*Received June 23, 1952. The research reported upon has been sponsored by the Geophysics Re- 
search Division of the Air Foree Cambridge Research Center under Contract AF 19(122)-217. 
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This may be regarded as a theoretical foundation of an empirical coagulation law valid 
for some cases of Brownian motion. Whytlaw-Grey (1932) and Sinclair (1950) discuss 
the range and validity of equation (6) in more detail. 

Of course, equation (6) holds exactly if ¢(z, y) = c. In this case, one can solve the 
main equation (3) under general initial conditions. Letting 2F (x, t) = cf(z, t), equation 
(3) becomes: 


IF (x, t) ie : ' -—" 
: = = | F(y, )F(x — y, t) dy — 2F(a, t) F(y, t) dy. (7) 


One introduces the Laplace Transforms: 
=~ /0 


Y(p, t) = | e "F(z, t) dz; Y(p, 0) = | e g(x) dx. (8) 


Then, using the standard properties of the Laplace Transforms, one obtains from equa- 


tion (7) the equation for y(p, @): 


,) f 
. EP: = v(p, t) — 2y(p, DWO, 0; ¥(p, 0) known. (9) 
Cl 
Its solution is: 
I (p, 0) : CN(0 
Vp, ) == >—e2” ve yoo=-—. 10) 
(Ni + 1) — a 0) 2 
Nie’! 
Now, inverting equation (8), one obtains the final solution: 
2 ] »y+ia pry / 0) 
(,) =a | =e — »; (11) 
: C(Nt + 1)° 2m J,-i 1 =" 0) 
vee’ 


y > 0 and sufficiently large. 
(From equation (11), one may obtain as a particular case a solution of equation (3) 
given by Schumann (1940), valid for a specific initial distribution only.) The conditions 
under which the above use of the Laplace Transform is justified are not discussed here, 
but it is easy to see from the physical situation that they are satisfied. 

A few particular results for initial distributions which may be of interest have been 


worked out. Two of them follow. 


If f(z, 0) = g(x) = A b(z — 1), 
then 
, 2N ~ Nt ) ,. CA 
(x = — > eee zr—(j + 1); N«a-—., 2 
We) = Cone + 1) §=0 (: + 1 lz — (9 DI; . 2 aad 


From equation (12) it follows, as is otherwise clear, that starting with single particles 
of equal size one will later have double, triple, m-tuple aggregates. The number of 
-~mp(t) 


m-tuples per unit volume at time ¢ is y(d)e 
If g(x) = f(z, 0) = A x" ee’, n integer, 


then 

F 2e° ** L Soi pi 

IM, ) = Nt + Dn + D aes exp (f La}, (13) 
where 


ae CAn! = enitA _ rite — 227 i/nt+1 
QN = aes b= [as | ; f=e . 


a 
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Equation (13) shows on closer analysis that as n increases the oscillatory terms become 
more and more dominant, giving progressively sharper and higher periodic peaks in 
F¢: t). 

Equation (3) was also solved under different assumptions for ¢(2, y) viz. ¢(z, y) = 
c(x + y) and (2, y) = cxy. The solutions, except for some special initial conditions, are 
involved, though quite straightforward. They will not be reproduced here, but may be 
obtained by transforming equation (3) by means of equations (4) and (5), introducing 
Laplace Transforms and solving the resulting partial differential equations. 

With physically more realistic assumptions for ¢(z, y), e.g. (2, y) = e(2'? + y'%)? 
(capture cross-section for two spheres), equation (3) becomes so involved that it has to 
be treated numerically. 

The following meteorological problem provides an illustration of the theory developed 
above. Marshall and Palmer (1948) found a simple expression for the distribution of the 
raindrop sizes. Their exponential relation, 

Ny = Nee*”, (14) 
where V,pdD is the number of drops in the diameter interval D — D + dD in unit 

value of space, 
N, is a universal constant, 
A is a rain-intensity parameter, 

fits the observed data of Laws and Parsons (1943) very well except at small diameters. 
This deviation at small diameters may be explained by assuming the distribution given 
by equation (14) subject to a collision process with the above form of ¢(z, y). The results 
are shown in figure 1. An attempt was also made to explain theoretically the whole dis- 


’ 
ic 


+ 














L 4 i 4 i 


s k 1 





me 
Fic. 1. The marked points represent the experimental data of Laws and Parsons, the straight line is 


obtained from the Marshall-Palmer relation, the curved line represents the calculated distribution. 
Np, A, D are explained under Eq. (14) in the text. 
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tribution given by (14), but the results are not yet satisfactory. This means probably 
that the form of ¢(x, y), chosen above, does not describe the physical situation suf- 
ficiently closely. 

The author is indebted to Professor J. S. Marshall, of McGill University, for sug- 
gesting the above problem, and to Dr. Walter Hitschfeld, of McGill University, for 


helpful suggestions concerning the physical model used. 
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VARIATION OF COEFFICIENTS OF SIMULTANEOUS LINEAR EQUATIONS’ 
By JOHN E. BROCK (Midwest Piping Supply Co., Inc.) 


1. Introduction. In dealing with sets of simultaneous linear equations it sometimes 
becomes necessary to modify one or more of the coefficients (to correct errors or for 
other reasons) after the solution has been obtained. Methods of obtaining the corre- 
sponding modifications in the solution with a minimum amount of computation have 
been treated by B. L. Weiner’, and in a discussion of Weiner’s paper (which bears the 
same title as the present paper), I. F. Morrison” indicates a matrix formulation of the 
analysis. In the present paper, we develop three particular prodedures each of which 
appears to be quite simple and quite useful. 

Matrices are used throughout in these developments. Both matrices and scalars will 
be denoted by ordinary italic letters. Letters with subscripts will denote matrix elements 
(scalars) and other scalars will be easily distinguished from matrices. A symbol denoting 
a matrix or matrix product when enclosed within parentheses and followed by subscripts 
= c. Super- 


will denote the appropriate element of the matrix, thus(ab);; = c;; where ab 
scripts will be used both to denote exponents and to serve as distinguishing indices; 
there should be no difficulty in distinguishing between these usages. 

It is presumed that we know the set of inverse coefficients for the system, that is, 
the matrix which is inverse to (reciprocal to) the original coefficient matrix. This is not 
a severe demand since the reciprocal matrix is frequently already known or can be com- 
puted without much trouble making use of calculations already performed in obtaining 
the original solution. 

So that the operations will not be obscured by voluminous calculations, the order, 
n, of the systems in the examples chosen for illustrative purposes is small. However, it 
is obvious that the advantage afforded by the procedures described here increases as n 
increases. 


1Received June 15, 1952. 
2B. L. Weiner, Variation of coefficients of simultaneous linear equations, Trans. ASCE, 113, 1349 








(1948) (Paper No. 2358). 
3]. F. Morrison, ibid., p. 1379. 
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2. Symmetrical Systems. Suppose a is a symmetrical, non-singular matrix and that 
we know a = a’, its inverse, and the solution, x, a column matrix, to the equation 


ar =. (1) 
Now let a be changed to 
a* =a+b, (2) 
where all elements of b are zero except only 
Dye = bee = B; (3) 


that is, the equal and symmetrically disposed elements a,, and a,, are increased by the 
amount 8. Our problem is to find a new solution 


marty (4) 
such that 
a*zx* =c. (5) 
In the derivation,‘ let the matrices b”* be defined by 
(b") is = 658; (6) 


where the (Kronecker) deltas are zero unless the upper and lower indices are the same 
in which case their value is unity. In establishing the relation (14) below, we use the 
“summation convention” which implies a summation from 1 to n inclusive over all re- 
peated indices except p and g which are fixed. Then 


(ab"),; = a, 55 , (7) 
(bab); = On, 8%5; , (8) 
b = Bb” + b”), (9) 
(ab);; = Bla@ipd; + ai955), (10) 
(bab) :5 = B (pd, 55 + agg 555; + Opp 555; + ary, 55 55), (11) 

(cebacd) 5; = B°(aipQanS; + ip qg 5; + Oighpp 5; + Oig&p_95) 
= Bayg(ab);; + B’(@ine aad; + aig@pp 5), (12) 

(babab);; = Batyg(bab);; + B°(aeerQtgqdt 5; + ApyQtqqd: 5; 


HF ApyQo¢ 5: 5; + AyqX%yp 5 5) 


I 


‘The reviewer has suggested adding at this point a remark concerning motivation. The general case, 
Section 4, came first to the writer’s attention. Success in summing the series (37) in particular cases leads 
to the analyses presented in Sections 2 and 3, Essentially one attempts to determine powers of the 
matrix m = ab, and in Section 2 is led to Equation (14) and in Section 3 to Equation (31); having these 
results, the summation is easily accomplished. However, in presenting the results in these two Sections 
(2 and 3), it is not necessary to exhibit the series. Once Equations (14) and (31) are established, the de- 
sired results follow simply by a little algebraic manipulation. Use of index notation, the Kronecker deltas, 
and the summation convention abbreviates the derivation of Equations (14) and (31). 
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Thus we have established the relation 


(Vol. XI, No. 2 


bmm = Abm + Bb, (14) 
where 
m=ab=a''b, (15) 
A = 28a,, (16) 
B = B(ay a, a.) (17) 
We now prove the important formulas 
m(mx) — (1 + A)(mz2) 18 
y = 
y i+i~e —_ 
and 
ne _ _ , m(mx) — (1 + A)(mz) 
a =; =; oa Q 
‘ee a 1+A-B (19) 
To do this, form the product a*z* and in simplifying note that 
am = aa b= b. (20) 
We have 
a*x* = (a+ b)(x+y) = art brt+ (act dy 
+l 
=e+ b+ j = : —% [mm — (1 + Aym)hz 
=c [1 A)b — Bb + bm bmm — (1 A)b — bm — Abm] ————— 
+ + + + bmm (1 + A)b br y. mT A-—B 
x 
=c+(l — Abm — Bb) —_———_;, = 2 
Cc (bmm {bm — Bb) ea es c (21) 


making use of Equation (14). Q. E. D. (Of course we require that B # 1 + A.) 


We may also easily obtain the corrected inverse a* 


in the form 


where 
mm — (1 + A)m 
+ ——_— tr, 
I+ i+ A-B8B 


= (a*)~'. Write Equation (19) 


If we postmultiply this matrix in succession by the columns of @ we obtain the corrected 


columns, i.e., the columns of a*. In other words 


a* = ya 
Example: 
i9 10 4] -1 2 0 
a=1}10 5 a}. a= 2 —3 —2]; z= 
4 2 | 0 -2 5 


2 ‘ 
—3|; c= |3). (24) 
wi 1 
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Now suppose 


” 10 4 ° 0 0 
a* = ; 10 5 Of. Then 6 =(-2):;0 O 1 and 
{ 0 1 , 1 O 
we calculate: 8B = — 2, p = 2,q = 3, A = (2) (— 2) (— 2) = 8, 
B = (—2)?{(—3)(5) — (—2)?] = —76; 1+ A—B = 85, 1+A=9, 
0 «i ] 
m=i0 } 6 mx = ;-—I18}; 
, —10 | el 
* . = | —140 
y= ima Us Aim | o = -5 ~~ 6] -18] + 85 =| 246] + 85 
ie 4 = 2 
¢ «=f «<§ 26 50 
"| 
a*=2r+y= —9; + 85. 
wi 


This is the easiest arrangement of the calculations. If the new inverse matrix is 
required, as would be the case, for example, if it were necessary to apply a second set 
of corrections, we would compute: 


85 40 20 


y =i © 5 —6] + 85; 
0 10 5 
—5 10 | "7 
a* = ya =| 10 —3  —40; + 85; z* =yrx =| —9; + 85. 
20 — 40 ; oa 


3. Single Corrections in Non-Symmetrical Systems. Suppose again that a, x, and C 
satisfy Equation (1), that a is modified as in Equation (2), and that we desire to find 
z* given in Equation (4) so as to satisfy Equation (5). However, now we do not require 


that a be symmetrical and we presume 6 has only one non-vanishing element 


B = 6,, . 
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Thus, using a notation introduced earlier, we can write 

















b = Bb", (26) 
m;; = (ab); = Ba:yp5; , (27) 
Mog = Baap 5 (28) 
Maedis = Baepd.5; , (29) 
(bm) ;; = Bard. 5° . (30) 
Thus, we have proved that 
bm,, = bm. (31) 
Now we argue that 
— MZ . 
i] = ser (32) 
y l + Mm. 
r* = ¥y2 (33) 
where 
m 
=][-—- = (34) 
4 1 + m,, 
for 
(a + b)m 
*¥7* = ar+ br - a 2 
a*a ax + 6a i+ x 2 
=c-+ [b+ bm,, — am — bm|-———— = ¢, (35) 
ie lL + Me 
making use of Equations (20) and (31). Also, as before 
a®* = ya. (36) 
Example: 
[1 41. 3] | —5 15 19 —8] T 1] | 2 
0-13 -!1 9 17 1 —12 —1 6 
a= a= +44; z= ; c=] 
3 10 2 1 10 —2 2 2 f 
11 —2 5 l | | 38 —31 —7 18 | 4 1 14) 











Suppose now that a;, = 3 is changed to 5. Then p = 3, g = 1, 8 = 2, 


| 19 0 0 0 


1 ue © © 
m=ab= + 22: Mog = 


19 
22 
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r22 0 0 0 
at «© € 
I - = = + 41 
. 2 0 41 O 
Ls £° 6-2 
l 22 | [ —5 s. 2 + =) 
—42 17 31 1 —22 
‘=7f = + 4l1; a* = ya = + 82 
84 7 20 —2 3 
| 48] Ll @ =f «7 31] 














4. General Case. Both the preceding procedures which give exact corrections in 
convenient finite form were originally motivated by the analysis which follows and 
evolved from it through success in actually summing the matrix series which appears 
below. 

Suppose again that a, x, and c satisfy Equation (1), that a is modified as in Equation 
(2), and that we desire to find x* given in Equation (4) so that Equation (5) is satisfied. 
We now make no restriction on b other than that the infinite matrix series 


2 3 9° 
y=l—-—-m+m-—-m-4+-:-:: (37) 
converge, and indeed in such a manner that its terms may be rearranged. 

Then, as before, Equations (22) and (24) hold, for 

a*t¥a* = (a+ byx =c+(-atayrt byt =, (38) 
since the quantity in parentheses may be written as 

at+ay + by) = -—a+a—am-+ am — am + am — 

+ b — bm + bm*® — bm* + 
= (b — am)\(I — m+ m’ — m*® + ---) = (b — amy, (39) 
by reference to Equation (20). 

In practice, to investigate the convergence and then to sum a sufficient number of 
terms of the series to assure the desired accuracy may prove prohibitively laborious. 
Therefore, it is suggested that the procedure implied in Equations (37), (22), and (24) 
be used only if it is clearly obvious that the corrections b are so small relative to a that 
not only is convergence assured but also adequate accuracy may be obtained from the 
first few terms of the series. 

Example: 


36 16 4 3 -8 6 17 o 


| 
-~ 
II 
| 
bo 
om 
| 
~~) 
bhom 
| 
bo 
wae 
8 
ll 
| 
or 
cm) 
| 
? 2) 
is) 
II 
| 
et) 
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0? 16 4 ? 0 0 
a® = | 15 9 2.99); b=;0 O —1] + 100; 
6 4 2 ’ 0 0 
3 0 4 002500 0 atl 
m=ab={!-6 0 —12]| + 1200 = | —.005000 0 —.010000/; 
3 0 —12 002500 0 pine 
21 0 | 000015 0 am 
m=)|]—54 0 —168; + (1200)° = | —.000037 0O —.000117); 
1 0 al 000031 0 al 
997515 O — .003291 
y=I—-m+m = .004963 I .009883 |; 
— 002469 0 pa, 


2.096271 


s*=yr= a 


7.049273 
The new inverse matrix a* may be calculated from Equation (24). A complete new 


solution gives the values 
2.096277 


z* = (ance 
pees 
and it is seen that good accuracy is obtained with our method using only the first three 


terms of the series given in Equation (37). 


ON A SOLUTION OF THE ENERGY EQUATION FOR A ROTATING 
PLATE STARTED IMPULSIVELY FROM REST* 
By RONALD F. PROBSTEIN (Princeton University) 

1. Introduction. The problem of the steady, laminar incompressible flow of a fluid 
over an infinite plate rotating at a constant velocity was first solved by von Karman 
[1]. Recently the associated heat transfer problem for the von Kaérm4én example was 

*Received Nov. 3, 1952. This work has been supported by the United States Air Force, Air Research 


and Development Command, under Contract AF 33(038)-250. 
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calculated by Millsaps and Pohlhausen [2]. Neglecting heat conductivity Nigam [3] has 
obtained the corresponding non-steady velocity distributions valid for small times for 
the rotating plate started impulsively from rest with a constant angular velocity. It is 
the purpose of the present paper to give the temperature distributions associated with 
the velocity fields found by Nigam. Apart from an academic viewpoint this solution 
offers an interesting mathematical analogy to the steady hypersonic laminar flow over 
a semi-infinite flat plate with a ‘‘self-induced’’’ pressure gradient. 

Under the assumption of an incompressible fluid a solution of the energy equation 
will be obtained for a thermally insulated surface and for a non-insulated surface with 
a wall temperature distribution which is a parabolic function of the radius. By analogy 
with Pohlhausen’s steady, laminar boundary layer analysis [5], the assumption of 
constant density requires that the azimuthal plate velocity must be small compared 
with the sound speed. Furthermore, in the non-insulated case as the viscous dissipation 
is to be included, it is implied that the energy dissipation associated with viscosity is 
of the same order as the energy associated with the imposed temperature difference. 
The maximum permissible difference between the wall and ambient temperature must 
therefore be small, that is, of the order of the square of the azimuthal plate velocity. 
Both of these conditions place a limit on the radius beyond which the solutions to be 
obtained are no longer valid. : 

Like Nigam’s solution of the equations of motion, the present solution of the energy 
equation will neglect all time dependent terms of order @ or higher; that is, only the 
early stages of the motion when ¢ is small will be considered. Such a procedure is tanta- 
mount to an asymptotic expansion in which the functions are expressed in ascending 
power series in ¢ and the first coefficients is obtained. This of course has the limitation 
that no information is given regarding the time to reach the steady state. (See Ref. 6 
for an analogous procedure in the case of the impulsively started cylinder.) 

2. Basic equations. The basic partial differential equations for the time dependent 
motion (@/dt ~ 0) of an incompressible fluid (density p = constant) are most suitably 
expressed for this problem in cylindrical coordinates r, 6, z with radial, azimuthal and 
axial velocity components w, v, w. Since the angular velocity of the plate 2 will be assumed 
constant for all time ¢ > 0 the problem is axially-symmetric and hence the partial de- 
rivatives with respect to @ vanish. Therefore the energy equation is written as 


{or aT aT. x (em 1 aT orl 
ote o 21 + + 


ot tu or 7? az | - p or” r or 
{_fau\? u\? dw/\’ av v\ (2) (2 aw) 
{9 “ ow = .* ov ou , gw 
"? (24) + of + (2 ) . (2 ) + dz * dz sa or / }’ (1) 


where, 7’ is the temperature of the fluid, c, the specific heat expressed in mechanical 
units, \ the thermal conductivity, and » the kinematic viscosity. 
Nigam has shown that if the infinite plate at time ¢ = 0 is suddenly made to rotate 


with a constant angular spin 2 about the axis r = 0, that the velocity components may 


1The term “self-induced” refers to the pressure gradient generated as a result of the curvature of 
the boundary layer over a flat plate; this pressure gradient in turn affects the boundary layer growth [4]. 
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be expressed as follows: 


u = Wrtf(n), (2) 
v = Qrg(n), (3) 
w= —4y'*0'*t"*h(n), (4) 


and the functions g(n) and f() satisfy the second order ordinary 


where, » = 2/2 (vt)? 
linear differential equations 


g’ + 2ng’ = 0, (5) 
f'’ + 2nf’ —4f = -—4¢°. (6) 


Note of course, that as the isothermal Blasius velocity distribution remained unaffected 
by heat transfer in the Pohlhausen analysis [5], so here, the solutions for the velocity 


components are unchanged by the heat transfer. 
It is now assumed that the energy equation is satisfied by a relation of the form 


e.T =c¢.T. + (Qr)*k(n) + Dvtm(n), (7) 


where 7... is the ambient temperature. On substituting Eq. (7) along with Eqs. (2)-(4) 
into the energy equation, neglecting all terms of order ¢ or higher, and equating coeffici- 
ents of ¢ and r’, the following two ordinary inhomogenous linear differential equations 


are obtained for the functions /(y) and m(n): 


| ed <i. 2nok’ = —ag’’, (8) 
m’’ + 2naom’ — 4om = —16k, (9) 
where ¢ = ¢,vp/d. The problem will now be particularized to the case of ¢ = 1 since 


an analytic rather than numerical solution is desired for the comparison of any results. 
The close relation between the above equations with ¢ = 1 and Eqs. (5) and (6) satisfied 
by g and f, is then immediately evident. 

3. Solutions of the equations. The solutions for m and k may now be expressed as 


m = 2(C + 2)g + (D— 1} — f], (10) 


k= 3[Cg+ D—(g—- 1), (11) 
where C and D are constants of integration, the additional two required integration 
constants A and B being contained in the general solution for f which is given by Nigam’s 
Eq. (10).* It is of course to be noted from Eq. (11) that the quantity involving & in the 
expression for the temperature distribution is a function only of the azimuthal velocity 
v. If g is required to satisfy the boundary conditions that either at time zero or infinitely 
far from the plate surface (7 = ©) the azimuthal velocity v = 0, and at the plate surface 
(n = 0) for all times greater than zero v = v,4;; = Qr, then g is given by 


g=1—27'” | e” dyn = 1 — erf n = erfe 7. (12) 
2The Prandtl number c,vp/d is equal to ac,/c, . In the case where cp = c, , the Prandtl number 


and o are of course identical 
The sign in front of the particular solution in Eq. (10) of Ref. 3 should be negative. 
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The particular solution of the energy equation for the thermally insulated plate is 
obtained by using the boundary conditions that the temperature gradient normal to 
the surface must vanish, and infinitely far from the plate surface and at time t = 0 the 
temperature has its ambient value 7. . The temperature distribution in this case is 
therefore 


¢,T = ¢,T. + v(wai1 — 30) + rtd | (xe *~ n erfe 7)” 
— 2n’ erfe » + (40)? ne} (13) 


For the case where the wall temperature is a parabolic function of the radius (i.e. 
Taii = constant r° + 7.) the boundary conditions for zero time and infinite distance 
from the plate are the same as in the insulated problem. On the other hand, the boundary 
conditions at the plate surface become m(0) = 0 and k(0) = a, where the dimensionless 
constant a controls the rate of change of the temperature gradient along the radius. 
Hence, the solution for the temperature is 


9,,\ ’ » : ° 2 
a) a sh + 214 {(e-e-” — 9 erfe 9)? 


tu 


eT =¢,T. +0 (! 


— (1 + 2a)n erfe n — w '7(1 + 2m’) erfe 9 
+ 9 ne" [Qn + (1 + 2a)]}. (14) 


The present approach can be extended to wall temperature distributions different 
from parabolic by use of a method similar to that adopted by Chapman and Rubesin 
[7] for solving the steady, compressible laminar flow over a semi-infinite flat plate with 
an arbitrary distribution of surface temperature. 

4. Discussion of results. The main point of interest in the solutions of the energy 
equation represented by Eqs. (13) and (14) is their mathematically analogous nature to 
the solutions for the flow over a semi-infinite flat plate with a self-induced pressure 
gradient. (See footnote 1 and Ref. 4). Examination shows that both of the solutions 
obtained can be considered to be partially composed of a Crocco type integral of the 
steady flow energy equation for the flow over a flat plate with a Prandtl number equal 
to unity and a constant surface temperature under the boundary layer approximations 
(i.e. T = a + bv + cv’). In addition there is what might be called an “induced time” 
effect represented by the Q’vt terms which correspond formally to the self-induced 
pressure gradient already mentioned. 

The portion of the solutions involving only the azimuthal velocity arises from the 
dissipation term »(dv/dz)* and the analogue of the constant surface temperature in the 
steady flow flat-plate case is the parabolic distribution chosen for the impulsively started 
rotating plate. 

By analogy to the work in Ref. 4 the terms in v and 2’vt are but the first of an asymp- 
totic expansion in ascending powers of ¢. Thus the present work for small times bears a 
closer analogy to the treatment in which the steady hypersonic flow over a semi-infinite 
flat plate for the region close to the leading edge is considered [8]. In this case solutions 
are obtained by expansion in powers of x! , rather than in powers of x as in Ref. 4 
(x is the distance along the plate surface with the leading edge as the origin). The latter 
solution, which is valid for the region toward the back end of the plate, would in effect 
correspond to the large-time asymptotic expansion for the present problem. 
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ON A CERTAIN ASYMPTOTIC EXPANSION 


By A. VAN WIJNGAARDEN (Mathematical Centre, Amsterdam) 


1. Introduction. This little note is concerned with the asymptotic expansion of a 
certain function, defined by a series resembling that of the modified Bessel function 


T(z), viz. 


Pp 
(2/2) 


(2) = ye (1,1) 


v 


This function plays a role in a certain physical problem, and values of /(z) for large 
positive values of the argument had to be determined. In itself the problem is of little 
general interest but it gives an opportunity to stress again the importance of the use of 
factorial series and to show the simple means they provide to derive an asymptotic 
expansion in such a case. 

2. The factorial series. The modified Bessel function /,(z) has the expansion 
(2/2)” 


(7) = “ ne 9 
1,2) ik + pr (2,1) 


k=O 


As k! k’’* ~ (k + 1/2)! this series for p = 1/2 is, apart from some trivial modifications 
as multiplication by (z/2)'’” and the addition of a term with k = 0, closely related to 
the series (1, 1). 


This relationship is used in the following way. The ratio (k + 1/2)!/(k! k'”*) can 
be expanded into a factorial series (cf. Nérlund, Differenzengleichungen): 


| ‘ ] *\p s+1/2) 

(—1) ‘ 

(k + 1/2)! 8 E 
— —— = - - ~ ‘ eo 
kik’? p> (k + 3/2)(k + 5/2) --- (k +8 + 1/2)’ > 0), 4,4) 


are generalised Bernoulli numbers. As these alternate in sign with in- 


where B,**'” 
creasing values of s, the terms of the series (2, 2) are all positive. Or, 


l s— 1/2 oe 1 
—5 = —)’ a aie ’ 5 ) (Bua 
kik” 2d ) ( 8 y (k +s + 1/2)! (k > 0), 3) 
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so that after inserting into (1, 1) and interchanging the order of the summations 


a r_\a{ 8 — l ?) (#+1/2) = (2/2)* 
fe = Lit ( . PP Deis +1! 


s=0 


_ 74? 1 Soe. dN 
ot ( 8 )e: ‘(3) oe eres? ox unh 


s=0 





or with regard to (2, 1): 


- : ’ a“ 1/2 (#+1/2) z tar = — i 
fe=> (-)'( ‘ Jp: ‘(2) T.41/2(2) oe ai: (2,4) 


s=0 


So far, no approximations have been introduced. It should be noted that the two 
terms between brackets in (2, 4) may not be separated as the latter term should give 
rise, according to (2, 3) to a divergent series. 

The functions /,,,,.(z) are, as well known, elementary functions. For the first few 


values of s they are respectively: 


I,,2(z) = (2/xz)'” sinh z, 


‘| cosh z — - ~ sinh :), 


I5/2(2) = (2/xz)' “la a” ; cosh z z+ :, sinh :), 


Zz 


~ 


15. 
Iz,.(2) = (2/xz)' — - © sinh z + — 2 cosh 2 — 73 sinh :), 


45. 105 105. 
*\sinh z — — a z + —sinhz — — cosh z + —j-sinh z}, 


1S 105 420 . 
I s3;9(2) = (2/x2) we z- sinh z + - cosh z — - sinh z 


€ 
x cosh z — ee sinh z). 


Moreover, the numerical values of the coefficients in (2, 4) between the sum-sign 
and the brackets are for the first few values of s respectively: 


8 rt = Elgon 


8 
0 ] 
I] 3/8 
2 65/128 
3 1225/1024 
4 131691 /32768 
5 4596669 /262144 


The series (2, 4) is convergent, and moreover for large values of z manageable. 








246 BOOK REVIEWS [Vol. XI, No. 2 


3. The asymptotic series. For practical purposes, it is however advantageous to 
give up the convergence of the expansion and turn it into an ordinary divergent asymp- 
totic series with simpler terms to be used for z > 1. To that end, one puts sinh z ~ cosh z 
~ e'/2, neglects the second term between the brackets in (2, 4), inserts the formulae 
given for J,,,,2(z) into (2, 4) and rearranges the terms with respect to decreasing powers 
of z. Then one obtains 


; -— :..,#.., 4, S71... , SR... ) 

(z) ~ e'(mz) ** +-2z ig Zz Zz ” eee & (3,1) 

J ( r 4% 139% + jog* T o048 8192 ‘a 
The given terms yield e.g. f(20) ~ 1.42507-10’; actually {(20) = 1.42508-10’. 


A note on my paper 


A GENERAL SOLUTION FOR THE RECTANGULAR AIRFOIL 
IN SUPERSONIC FLOW 
Quarterly of Applied Mathematics, XI, 1-8 (1953) 
By JOHN W. MILES (University of California) 
It appears that the application of the Lorentz transformation to unsteady flow prob- 
lems was suggested previously by P. A. Lagerstrom (unpublished lectures at California 
Institute of Technology, Pasadena, 1948) and by C. Gardner in his paper 7'%me-dependant 


linearized supersonic flow past planar wings, Comm. on Pure and Appl. Math. 3, 33-38 
(1950). 


BOOK REVIEWS 


Electronic analog computers. By G. A. Korn and T. M. Korn, McGraw-Hill Book 
Company, Inc., New York, Toronto, London, xv + 378 pp. $7.00. 


This book is both 1 text and an engineering reterence book on the subject of d.e. electronic analog 
computors. The type of computer under consideration is that employed by many commercial and govern- 
ment institutions. These computers are used in the solution of such problems as are encountered in 


the design of feedback control systems and in the development of new aircraft configurations. 


In the book, the problems associated with the use of computers are first discussed. Examples are 


given of setup prt cedures and applications to show both the limitations and the versatility of these 
devices. Following this, the principal engineering design considerations associated with the computer 


components are discussed. The d.c. amplifiers employed are analyzed at length as are the various means 
of obtaining special functional representations. The book is concluded with a description of several 
of the commercially available computers. 

As a text on electronic analog computers, it is unfortunate that the a.c. computers are not considered 
as well as the d.c. computers. Because of this the book is not sufficiently balanced for a text on analog 
computers and would have to be supplemented by additional reading. The principal value of the book 
is for the orientation of mathematicians and others concerned with the solution of dynamic problems 
on analog computers that have either been specially designed or commercially procured. 


H. T. Marcy 




















THE FOURTH SYMPOSIUM ON PLASTICITY 


will be held at Brown University, Providence, Rhode Island on 
September 1, 2 and 3, 1953. The three principal topics of this 
Symposium will be: 


Stress-strain relations (perfectly plastic, work- 
hardening, visco-elastic, and visco-plastic materials) 
Distribution of stress and strain in bodies exhibiting 
inelastic properties (static and dynamic problems, 
infinitesimal and finite deformations) 


Structural analysis and design in the inelastic range 
(limit analysis, static and dynamic loading, shke-a 
down, inelastic buckling). 


Both foreign and American experts in the field will present papers 
covering their most recent experimental and analytical research, 
and everyone is invited to participate in the discussion. 

Certain sessions of the Symposium, containing papers of par- 
ticular interest to civil engineers, are being sponsored, jointly with 
Brown University, by the Engineering Mechanics Division of 
the American Society of Civil Engineers. The remaining sessions 
are under the joint sponsorship of the Office of Naval Research and 
Brown University. 

Persons interested in receiving further information concerning 
the Symposium should write to Professor H. J. Weiss, Graduate 
Division of Applied Mathematics, Brown University, Providence 12, 
R. I. Copies of the program and reservation cards will be sent out 
approximately six weeks before the date of the Symposium. 
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Introduction to the Foundations of Mathematics 

By RAYMOND L. WILDER, University of Michigan. A critical inquiry into the nature of 
mathematics and the foundations upon which it is built. This unusual book clearly explains: 
the source and evolution of mathematics . . . what constitutes mathematics and mathe- 
matical existence . . . contradictions and paradoxes and methods of avoiding them. . . the 
origin and growth of present day theories. 

‘The book is filled with a vast amount of information organized so that it can be under- 
stood. Expounding some of the most difficult material in the world, it reads smoothly. I 
found that almost all of the questions I was provoked to ask, or objections that occurred to 


me, were attended to by the book in due time.’’—Dr. Leo Zippin, Queens College. 
1952. 305 pages. $5.75. 


Numerical Solution of Differential Equations 
By WILLIAM EDMUND MILNE, Oregon State College. Numerical methods for solving 
partial and ordinary differential equations are treated thoroughly for the first time in 
English in this combination computer's manual and explanatory treatise. Every important 
topic is illustrated with specific numerical examples, demonstrating exactly how to apply 
processes to actual problem situations. 
“*May I express my deep appreciation . . . for this most useful volume which will be of great 
value to me in my professional work.’’—Herbert E. Salzer, National Bureau of Standards. 
A publication in the Wiley Applied Mathematics Series. I. S. Sokolnikoff, Editor. 

1953. 275 pages. $6.50. 


Demand Analysis: A Study in Econometrics * 


By HERMAN WOLD, University of Uppsala, in association with LARS JUREEN, Stare 
Agricultural Marketing Board of Sweden. A systematic account of demand analysis methods, 
employing for illustrative material the empirical studies of the authors. The book provides 
a clear exposition of the traditional method of least squares regression, presents new lines 


of theoretical argument, and many new types of problems. 
1953. 358 pages. $7.00. 


Stochastic Processes * 
By J. L. DOOB, University of Illinois. The mathematical background to the advanced theory 


of probability, including material never before published in English. 
1953. 654 pages. $10.00. 


50-100 Binomial Tables* 
By HARRY G. ROMIG, Hughes Aircraft Company. This important new work provides tables 
for the positive binomial (¢ + p)", where q = 1 — p. It covers the range of m values from 


50 to 100 in steps of 5 and the range of p values from .01 to .99 in steps of .01. 
1953. 172 pages. $4.00. 


An Introduction to Linear Programming 
By A. CHARNES, W. W. COOPER, and A. HENDERSON, all at the Carnegie Institute of 
Technology. A valuable guide to the theory and mathematics of a rapidly expanding field in 


industrial efficiency and methodology. 
1953. 74 pages. $2.50. 


* Wiley Publications in Statistics. Walter A. Shewhart, Editor. 


Order on-approval copies today 


JOHN WILEY & SONS, Inc. 440-4th Ave., New York 16, N. Y. 


























